2015年全国大学生数学建模竞赛国家一等奖论文-A题-太阳影子定位模型的分析

word文档可自由复制编辑word文档可自由复制编辑太阳影子定位模型摘要趋势；对杆长呈正比关系增长；对经度呈现先急剧增长到峰值再突变为01中提供的数据，画出了天安门广场上直杆的太阳影子分布曲线图。20 21? i i

? 归i 归i|DA-DA'|,min |S

-S'|。i=1

i=1)杆长为，太阳方向角残差比为1.8%，影长残差比为型求解的时间复杂度较高。为提高计算速度，引入了粒子群算法。分别对附件2和3中的数据进行分析，确定出的地点坐标分别为)，)和))2为西藏阿里，日期为8/14或4/，附件3为广西梧州市，日期为12/27或12/14。可以从视频中直接得到的影长并不是实际长度，而是其投影长度，这里采用基于Hough变换和透视变换的图像矫正法，对斜视图像进行矫正，得出实际影长。然后将得到的数据带入问题二的模型中，给出视频拍摄地点为)题三的基于粒子群算法的目标规划模型，求解出视频拍摄地点为)——内蒙包头市境内，拍摄日期为6/11或7/13。建筑物群合理间距问题。关键词：双目标规划粒子群算法Hough变换透视变换一、 问题的重述1.1问题的背景与太阳活动有着密切的联系。而古代智慧的先民就利用了这个现象制作了日晷，时间和地理信息，是非常方便可靠的。1.2问题的提出的地点和日期的一种方法。针对上述背景和应用需求，提出以下问题：1.建立影子长度变化的数学模型，分析影子长度关于各个参数的变化规律，并应用你们建立的模型画出2015年10月22日北京时间9:00-15:00之间天安39度54分26秒,东经116度23分29米高的直杆的太阳影子长度的变化曲线。1个可能的地点。3.根据某固定直杆在水平地面上的太阳影子顶点坐标数据，建立数学模型2和附件3的影子顶点坐标数据，给出若干个可能的地点与日期。4．如果已有一根直杆在太阳下的影子变化的视频，并且已通过某种方式估计出直杆的高度为2与日期。二、 问题的分析2.1问题一的分析过程。2.2问题二的分析值之和达到最小时，所得经纬度即为杆子的地点坐标。2.3问题三的分析与问题二不同的是，该问中日期是个未知量。首先考虑沿用上一问的模型。365理的适应度函数，期望提高计算效率。2.4问题四的分析得出实际影长。然后再分别利用问题二和问题三的模型进行求解。三、 问题的假设1.2.地球是一个完美的球形，不考虑海拔、地球扁率的影响。3.无光线衍射造成的影子减淡现象。5.在小尺度考虑直杆投影问题时，地表是绝对水平的。6.不考虑地球公转的影响。7.题目所给的数据是真实的，可靠的。四、 符号说明符号说明q太阳高度角L水平地面上直杆长度S水平地面上直杆影子长度h当地地方时时角d当地太阳赤纬当地纬度值g当地经度值t当地地方时t0北京时间n当日日期序号P当地经纬度坐标A当地太阳高度角T时间（包含时刻与日期）b离散化后的方位朝向ei每组数据影长与x坐标轴的夹角xi，yi题目附件提供的第i组直杆坐标d日期注：其他符号将在下文中给出具体说明。五、 模型的建立与求解5.1问题一5.1.1模型的建立一些地理学定义，以及重新定义一些本模型需要用到的参数。上，直杆长度与影长的比值即为太阳高度角的的正切值： tanq=L S图1太阳方位角示意图其中L为杆长，S为影长。通过查阅，得知太阳高度角的计算公式为： sinq=coshcosdcosf+sindsinf (2)其中q为太阳高度角，h为地方时时角，d为当时的太阳赤纬，f为当地纬度。以一个地方太阳升到最高的地方的时间为正午1212时之间等分为24个小时，所成的时间系统，称为地方时。地球上每一个地点都有角时，要先将北京时间换算为当地地方时t：g-120?t=t0+

()其中t0为北京时间，g为当地经度。根据某地地方时，可以换算出当地的地方时时角。地方时时角h即为当地与子午线之间相差的角度：h=15按(t-12) (4)纬的最高度数为'，夏至时太阳的赤纬为，冬至时太阳的赤纬为-23?26'。春分和秋分时太阳的赤纬为0?。d可用下面这个公式来计算：d=23.45?sin?(284+n)? (5)? 365 ?? ?其中n为当日日期序号，1月1日时n=1，以此类推得10月22日n=联立式子(2)-(5)得到方程组：? Ltanq=S?S??sinq=coshcosdcosf+sindsinf?t=t0+?

g-120?15?

(6)?h=15按(t-12)? ? + ?=23.45?sin?(284 n)?求解上述方程组，得：

? 365 ?S=L/tan(arcsin(cos(15按(t-12))cos(d)cosf+sin(d)sinf)) (7)可见，影子长度的变化与当地地理位置P(j,f)，直杆长度L，时间T(t,d)这三个参数有关。5.1.2模型的求解首先计算题目所给条件下的h与f(2)9:00-15:00制成下表以供参考：表1影子长度分布北京时间太阳高度角/度影长/m9:0021.187.7412:0037.883.8512:1437.993.8415:0025.356.33将影长随时间变化的情况用绘制成图像：图2影子长度变化曲线2中已经给出；对于地点，将其分为经度g和纬度f这两个参数进行运算求解，求解后的变化规律图如下：图3影长与日期关系曲线图4影长与杆长关系曲线图5影长与经度关系曲线图6影长与纬度关系曲线5.1.3结果分析从图3最小点出现在12:14附近时，会产生时差，显然是符合常理的。从图4tanq=LS了模型的准确性。从图5地平线以下，就不存在影子，影子长度也为0。从图65纬度的增加太阳高度角减小，从而导致了影长的增加。5.2问题二——基于网格逼近算法的双目标规划模型5.2.1数据预处理x轴的方位朝向角为b(0??b<360?)轴的方位朝向角为0?。21取部分呈现在下表中：表3影子方位角数据表北京时间x坐标/米y坐标/米方位角/度14:421.03560.4973b+25.63?14:571.20870.5255b+23.49?15:121.39550.5541b+21.65?15:271.60330.5833b+20.02?15:421.82770.6135b+18.55?5.2.2模型的建立定义为从北方沿着地平线顺时针量度的角：cosA=sind?cosf-cosh?cosd?sinfcosa

(8)上述公式可以用来计算近似的太阳方位角，不过因为公式是使用余弦函数，180180(下午)，方位角应该大于180度，即要取补角的值，故作如下修正：?arccos?sind?cosf-cosh?cosd?sinf?,

h<0? ? cosa ??A=? ? ???360?-arccos?sind?cosf-cosh?cosd?sinf?,h?0

(9)? ? cosa ?? ? ?由上述推导可得，方位角A与坐标系方位朝向b的关系为：=b+ei,

i=,21 (10)其中ei为每组数据影长与x坐标轴的夹角：ie=arctanyi,i

i=2, ,21 将上一问题的某些参数关系与式-列成方程组：??ei?

=arctanyi,

i=2, ,21?=b+ei,?

i=2, ,21? + ?? ?=23.45? ?(284n)365 ?

12)? ? ??h=15按(t-12)?=cos(15按(t-12))cos(d)cosf+sin(d)sinf?t=t0t=t0+? 15?将上述求太阳高度角的方程化为如下关系式：cos(A)=

(sin(d)?cosf-cos?按(t+15 15 0

g-120?15?

-12)??cos(d)?sinf)

? ? g-120? ??cos(arcsin(cos?15按?t0+ -12??cos(d)cosf+sin(d)sinf))? ? ??下面对各个参量进行影响因素的分析：g-120?t=t0+

，h=15按(t-12)?(284+n)?当地经度g的变化以及时间t0有关；对于赤纬角，由d=23.45?sin? ?? 365 ?为赤纬角dsinq=coshcosdcosf+sindsinf式子可关，我们可以将其理解为地点的变化；对于太阳方位角，由式子(8)可知，太阳大小有关，同样可以将其理解为地点参数的变化。综上，我们对参量的影响因素进行总结，可以得到下表：表4影响因素分析参量影响因素时角(h)经度g，时间t0赤纬角)常量太阳高度角)地点坐标(g,f)，时间t0太阳方位角(A)地点坐标(g,f)，时间t0可以发现，在式中，方程右端的未知量只有地点坐标(g,f)和北京时间t0，故可以将其简化为如下形式：F,f,t0)=(14)式中，映射F表示方程右端关于地点坐标(g,f)这两个自变量的函数。word文档可自由复制编辑21个不同时刻对应的太阳方位角(i=2,...,21个不同的方程，将其列成方程组如下所示：?F,f,t0)=cos()??F,f,t0)=cos()??,f,t0)=cos()由式(9)(g,f)计算得出此处不同时刻的太阳方位角(7)' 'i据地点坐标(g,f)算出此处不同时刻下的影子长度S'，对其进行归一化处理，可i以得到：

'Si-'

1211'?Si'S='S=归i S'

21i=1S'

16)imax

imini对此我们可采用遍历算法对不同经纬度(g,f)下的太阳方位角的插值DA'以i及归一化后的影子长度S归i优解。

进行求解，从而与实际值作差进行比较，进而得出最故此，我们建立多目标规划模型：?20?? ? i i?目标函数：??|-|i=121

17)?min |S

-S'|? ? i=1?-180??g?180??约束条件：?-90??f?90????0??q?90??

(18)式中，表示实际的太阳方位角，S归i表示归一化后的实际影长。5.2.1模型的求解0.01，对其进行编程求解：表5定点坐标数据经度纬度杆高19.12?表6误差分析表太阳方向角残差影长残差太阳方向角残差比影长残差比0.0260.014word文档可自由复制编辑即杆子所处地点坐标为)，在海南东方市境内。5.2.1结果分析)的影子长度残差比为，太阳方向角残差比为1.8%明得到的经纬度坐标)是较为精确的。5.3问题三5.3.1模型的建立变的情况下，要求模型输出杆子所处的地点和日期。由表3可以看出，日期影响的参量为赤纬角，其他参量并不受日期的影响，太阳赤纬公式为：d=23.45?sin?(284+n)? (19)? 365 ?? ?式中，n为当日日期序号，1月1日时n=1，以此类推，日期每增加一天，对n进行加1即可。n及经纬度坐标(g,f)进行遍历求解，建立与上一问相似的多目标规划模型，从而求出最优解：?20?? ? i i?目标函数：??|-|i=121

(0)?min |S

-S'|? ? i=1?-180??g?180?约束条件：??-90??f?约束条件：???0??q?90??>0

(21)式中，表示实际的太阳方位角，S归i表示归一化后的实际影长。5.3.2模型的求解编程求解：表7杆子所处地点和日期及误差分析数据表附件二附件三经度经度纬度纬度杆高杆高太阳方向角残差比太阳方向角残差比影长残差比影长残差比日期8/14或4/日期/27或12/145.3.3模型的改进——基于粒子群算法的目标规划模型对目标函数进行求解，从而降低模型求解的时间复杂度。假定有一个D维的目标搜索空间有n个微粒组成了一个粒子群其中每个微粒都用一个D维的向量描述，将它的空间位置表示为i=(1,i2,,D),i=,,,n这可看做目标优化问题中的一个解代入适应度函数计算出适应度值可以衡量微粒的优劣；第i个微粒的飞行速度也是一个D维的向量，记为i=(1,i2,,D)；第i个微粒所经历过的具有最好适应值的位置称为个体历史最好位置记为i=(1,i2,,D)整个微粒所经历过的最好位置称为全局历史最好位置，记为g=(g1,g2,,D)，粒子群的进化方程可描述为(t+1)=(t)+(t)((t)-(t))+(t)(pgj(t)-(t))(t+1)=(t)+(t+1)

(2)其中i表示第i个微粒，j表示微粒的第j个维度，t表示第t代，c1，c2个加速常量通常的取值范围是(,2)，1的随机函。

U(0,1)

U(0,1)是两个相互独立可以调节微粒去往自身周边的最好位置，c2可以调节微粒去往整个粒子群所能找到的最好位置。散的坐标值点：mi=(g,f,d)“粒子Zj了能够寻找到足够好的优化解，而且避免计算量的增大，本模型中设置了20个“粒子”来进行优化解搜索，即令i=。由于本问中直杆坐标和日期都是未知的，因此将粒子设为储存三个参数：?Z1=经度g???Z2=纬度f (23)??Z3=日期dg=rand(-180,+180)，f=rand(-90,+90)，d=rand“最优解”很可能会是局部最优解，而非全局最优解。度值Q做定义为：? ?1?2?R=? ??Q=? ?V??

(24)=

20 21|DA-DA'|+ |S

-S'|?? i=1

i i ? 归i 归ii=1记录做比较。如果能够比最好记录更好就将此值设为最好记录。接下来，参照式(22)，对每个粒子的去向进行调整。每次调整时，都要对?(0,2)，?(0,2)，?(0,1)，?(0,1)进行重新随机取值，保证粒子行为的随机性：既能有向最优解靠近的趋势，也能存在寻找其他更优解的可能。时间寻找解。这里令迭代次数为1000。5.3.4粒子群算法模型的求解表8杆子所处地点和日期及误差分析数据表附件二附件三经度经度纬度纬度杆高杆高太阳方向角残差比太阳方向角残差比影长残差比影长残差比日期8/14或4/日期/27或12/145.3.5模型的综合结果分析件12到的日期分别关于夏至日6/22及冬至日12/对称，这是由于太阳的运动是关知参量增加时，可以使用粒子群算法进行求解。5.4问题四5.4.1数据的采集及预处理——基于Hough变换及透视变换的图像矫正过的相关知识从图片中还原出影子的真实长度及变换的角度。透视投影的一个主要特征是可以延伸到无穷远的物体在照片平面中可以获点方向。而消隐点在水平方向上的延长即为地平线。据此，可以引出地平线的性质：在透视图里，从地平线上任一点引出的两条直线都是平行的。到消隐点，消隐点在水平方向的延长线就是地平线。令地平线与直杆的交点为P则P点就是直杆向其正前方的消隐点从直杆底部引出一水平辅助线B,则从P点到直杆底部O的连线L与水平线垂直。将影子端点R与PAP点引出的A平行于LA与B也垂直。因此，?RQO为直角，三角形DRQO为直角三角形。下面利用Hough变换[5]提取平行直线，我们知道，当图像中存在很多条直线若直线的宽度大于1，则舍弃。这个过程可以表示成：

m=ax+by+c,ox+py+1

(5)n=dx+ey+f (5)式中，斜视图像点的像素坐标为(x,y)，正视图像的像素坐标为(m,n)，a,b,c,d,e,f,o,p,q为透视变换参数。先，利用对其进行图像提取（提取速度为每2产生的22幅图片进行灰度处理和高斯滤波处理:灰度处理高斯滤波中阴影也变得更加清晰。5.4.2模型的建立1.在拍摄时间已知的情况下：代入到问题二的多目标规划模型中，即可得到视频拍摄地点(g,f)。?20?? ? i i? 目标函数：??|-|i=121

(6)?min |S

-S'|? ? i=1?-180??g?180?? 约束条件：?-90??f?90????0??q?90??

(27)式中，表示实际的太阳方位角，S归i表示归一化后的实际影长。2.在拍摄时间未知的情况下：间复杂度，并且误差比很小。题三的目标规划模型中利用粒子群算法进行求解：?20?? ? i i? 目标函数：??|-|i=121

(8)?min |S

-S'|? ? i=1?-180??g?180? 约束条件：??-90??f? 约束条件：???0??q?90??>0

(29)式中，表示实际的太阳方位角，S归i表示归一化后的实际影长。5.4.3模型的求解表9杆子所处地点和日期及误差分析数据表已知日期未知日期经度经度纬度纬度杆高杆高影长残差比影长残差比日期——日期6/11或7/135.4.4结果分析地点坐标为))经6频日期为7/13，两者刚好关于夏至日6/22对称，验证了模型的可靠性。六、 模型的检验6.1误差分析1.标准解与模型所求解的误差：由于本模型使用了目标规划方法对问题进发现其最大值为2.局部最优解造成的地理坐标误差：粒子群算法的本身原理会导致其收敛过程中，有时候需要人为干预来筛去一些的地区，可能造成冗余解或缺解。3.视频数据预处理过程中的误差：利用Hough变换寻找出的平行线并不一定是处在同一水平面上的，因此做透视变换的时候可能会产生误差。6.2敏感性分析我们是以每2分钟一帧的速度，从0时刻开始对其提取了22组数据。如果将起始时刻进行变换，即可以取1,5,...,39这些时间节点的视频数据，在已知日期的情况下带入模型进行求解，从而检验视频数据的改变，对模型结果产生的影响，得到的结果如下：表10敏感性检验对比数据原始数据检验数据经度经度纬度纬度健性。七、 模型推广可以应用到求建筑物群合理间距问题。平行的铅垂面,其在地面上的投影为两平行线,当所取的轨迹线上的时间差刚好等于要求的必要日照时间数时,两幢建筑之间的垂直距离即为合理的间距,这种方法适用于两幢建筑物等高的情况。假设楼体1的投影直线方程为=2的投影直线方程为y2=kx2+b，需要获取的日照时间数为ts1与投影曲线的交点(x0,y0)? 2+2?型可求解太阳高度角q=arccot??y0L

?h?? ?t=h +12(xt=h s 15?楼间距d=

b =1+k2

y'-kx'。1+k1+k2模型的优点：1.中需要对模型代码进行部分修改就可以进行计算。2.搜索能力，一般不会陷入局部最优。3.实用性高。只需要一段描述日影变化过程的视频以及相应的时间就可以求出当地的经纬度坐标。模型的缺点：1.只能应用在水平地面上。没有考虑斜坡和坎坷地面上的情况。2.对于拍摄角度不好的视频，该模型产生的误差较大。3.对形状规格较为复杂的建筑，本模型不能适用。九、 参考文献[1][2]袁新生，邵大宏，LINGO和EXCEL在数学建模中的应用，北京：科学出版社，2007。[3]谢金星，薛毅，优化建模与LINDO/LINGO软件，北京：清华大学出版社，[4]程序设计与应用，高等教育出版社，2012。[5]孙荣峰，快速霍夫变换算法，计算机学报，24(10)：1104-1109，2001。[6]郑鹏飞，林大钧，基于影子轨线反求采光效果的技术研究，华东理工大学学报(自然科学版)，36(3)：458-463，2010。[7]何援军，透视和透视投影变换，计算机辅助设计和图形学学报，17(4)：734-739，2005。[8]武琳，基于太阳阴影轨迹的经纬度估计技术研究，天津大学，2010。[9]叶德意，基于自适应变异粒子群算法的分布式电源选址与容量确定，电网技术，35(6)：155-157。十、 附录10.1数据表格表11太阳高度角与影长数据时间太阳高度角影长921.187959167.7392853859.122.085392227.3935161449.222.966792967.0789613349.323.831452266.7918104359.424.678639376.5288782479.525.507602246.2874837049.626.317568086.0653559259.727.107744185.8605606129.827.877318955.6714418429.928.625463295.4965756091029.351332225.33473245310.130.05406695230.732796985.04599397110.331.38664334.91736631310.432.014721014.79825987810.532.616143074.68805864210.633.190024174.58622312410.733.7354854.49228050310.834.251656914.40581626710.934.737686894.326467111135.192742924.25391487811.135.616019514236.006743424.12812394911.336.364179584.07443159911.436.687636934.02662174911.536.97647433.98453738511.637.230106033.94804464211.737.448007373.91703073811.837.629719483.89140223811.937.774853963.8710836071237.883096753.8560160412.137.954211353.84615654912.237.988041333.8414772912.337.9845123.841965126word文档可自由复制编辑12.437.943631143.84762142512.537.8654893.85846207812.637.750257273.8745177512.737.598187323.89583435812.837.409607513.92247378912.937.184919883.9545148611336.924596093.99205454613.136.629172764.03520946313.236.299246474.08411768413.335.935468274435.538537974.19986668813.535.109198314.26711188913.634.648229054.34092552713.734.156441154.42159297613.833.6346714.50944049913.933.083774964.6048405811432.504624074.70821816914.131.898099134.82005799414.231.265086084.94091318614.330.606471845.07141548114.429.923140445.21228735114.529.215969595.36435651114.628.485827685.52857336414.727.73357115.70603210914.826.960041975.8979964414.926.166066196.1059310611525.352451896.331540604表12视频采集数据时间节点影长12.42402402422.39108108132.35783783842.32483483552.31582582662.27654654772.24654654782.20741741792.176066066102.162462462112.121861862word文档可自由复制编辑122.105255255132.09039039142.054144144152.027087087162.003063063171.966996997181.94191.912942943201.900900901211.876876877221.85885885910.1源程序第一问程序：clc;clear;shadow1=zeros(1,601);%初始化不同时间段影长sthet=zeros(1,601);%初始化不同时间段太阳高度角的sin值thet=zeros(1,601);%初始化不同时间段太阳高度角n=295;%1月1日到10月22日间的天数time=[9:0.01:15];%从9点运算到15点，步长为0.01hh=[-48.60833333:0.15:41.39166667];%9点到15点地方恒星时系统下的时角delta=23.45*sin(2*pi*(284+n)/365);%太阳赤纬height=3;%旗杆高度phi=39.9072222222222;%天安门纬度fori=1:601%从9点运算到15点，共运算601个值sthet(i)=cos(h(i)*pi/180)*cos(delta*pi/180)*cos(phi*pi/180)+sin(delta*pi/180)*sin(phi*pi/180);%计算出太阳高度角的sin值thet(i)=asin(sthet(i))*180/pi;%计算出太阳高度角shadow1(i)=sqrt((height/sthet(i))^2-(height)^2);%计算影长endplot(time,shadow1)%输出图形ylabel('影长/m');xlabel('时间/时');title('影子长度与时间的关系');gridon;第二问程序clcclearall%附件一数据x=[1.03651.06991.10381.13831.17321.20871.2448...1.28151.31891.35681.39551.43491.47511.516...1.55771.60031.64381.68821.73371.78011.8277];y=[0.49730.50290.50850.51420.51980.52550.5311...0.53680.54260.54830.55410.55980.56570.5715...0.57740.58330.58920.59520.60130.60740.6135];length=sqrt(x.^2+y.^2);sOLen=sum(length);Length=length;length=(length(2:21)-length(1:20))./(length(2)-length(1));angle=acos((x(2:21).*x(1:20)+y(2:21).*y(1:20))./(sqrt(x(2:21).^2+y(2:21).^2).*sqrt(x(1:20).^2+y(1:20).^2)));result=[]; %经纬度t=14.7:0.05:15.7; %北京时间min=0.05;kkk=0.1;lon=80:kkk:120;beginl=10*ones(1,size(lon,2));endl=40*ones(1,size(lon,2));n=108;%2015年4月18日delta=23.45*sin(2*pi*(284+n)/365); %赤纬角fori=1:size(lon,2)ts=t-(120-lon(1,i))/15;%当地地方时omega=15*(ts-12); %时角beginl(1,i);endl(1,i);forlat=beginl(1,i):kkk:endl(1,i)phi=lat; %纬度h=asin(sin(phi*pi/180)*sin(delta*pi/180)+cos(phi*pi/180)*cos(delta*pi/180)*cos(omega*pi/180));%太阳高度角 a=2*pi-acos((sin(delta*pi/180)-sin(h)*sin(phi*pi/180))./(cos(h)*cos(phi*pi/180)));%太阳方位角l=100./tan(h);%影长l=(l(2:21)-l(1:20))./(l(2)-l(1));a=a(2:21)-a(1:20);% if(sum(abs(l-length))<min)% ganH=mean(Length.*tan(h));% result=[result;lon(1,i)latganHsum(abs(l-length)) n];if(sum(abs(l-length))<min)%寻找最小点min=sum(abs(l-length)); %更新min为两者更大的ganH=mean(Length.*tan(h));result=[lon(1,i)latganHabs(l-length)abs(a-angle)n];%最小点保存在result内endclearh;endn

endiplot(result(:,1),result(:,2),'r*')function[ParSwarm,OptSwarm]=BaseStepPso(ParSwarm,OptSwarm,ParticleScope,MaxW,MinW,LoopCount,CurCount)ifnargin~=7 %输入容错error('error')endifnargout~=2 %输出容错error('error')end%线形递减策略w=MaxW-CurCount*((MaxW-MinW)/LoopCount);%得到粒子群群体大小以及一个粒子维数的信息[ParRow,ParCol]=size(ParSwarm);%得到粒子的维数ParCol=(ParCol-1)/2;SubTract1=OptSwarm(1:ParRow,:)-ParSwarm(:,1:ParCol);%求解出历史最优值与当前位置的差值c1=2;c2=2;%完成一次粒子位置速度最优值的更新迭代forrow=1:ParRowSubTract2=OptSwarm(ParRow+1,:)-ParSwarm(row,1:ParCol);%计算出全局最优值与当前该粒子位置的差值%速度更新公式TempV=w.*ParSwarm(row,ParCol+1:2*ParCol)+c1*unifrnd(0,1).*SubTract1(row,:)+c2*unifrnd(0,1).*SubTract2;%限制速度的代码forh=1:ParColifTempV(:,h)>ParticleScope(h,2)TempV(:,h)=ParticleScope(h,2);endifTempV(:,h)<-ParticleScope(h,2)TempV(:,h)=-ParticleScope(h,2)+1e-10;%加1e-10防止适应度函数被零除end

end%更新该粒子速度值ParSwarm(row,ParCol+1:2*ParCol)=TempV;a=0.729;%约束因子TempPos=ParSwarm(row,1:ParCol)+a*TempV;%限制位置范围的代码forh=1:ParColifTempPos(:,h)>ParticleScope(h,2)TempPos(:,h)=ParticleScope(h,2);被零除

endifTempPos(:,h)<=ParticleScope(h,1)TempPos(:,h)=ParticleScope(h,1)+1e-10;%加1e-10防止适应度函数end

end%更新该粒子位置值ParSwarm(row,1:ParCol)=TempPos;%计算每个粒子的新的适应度值ParSwarm(row,2*ParCol+1)=AdaptFunc(ParSwarm(row,1:ParCol));ifParSwarm(row,2*ParCol+1)>AdaptFunc(OptSwarm(row,1:ParCol))OptSwarm(row,1:ParCol)=ParSwarm(row,1:ParCol);end

endword文档可自由复制编辑%for循环结束%寻找适应度函数值最大的解在矩阵中的位置(行数)，进行全局最优值的改变[maxValue,row]=max(ParSwarm(:,2*ParCol+1));ifAdaptFunc(ParSwarm(row,1:ParCol))>AdaptFunc(OptSwarm(ParRow+1,:))OptSwarm(ParRow+1,:)=ParSwarm(row,1:ParCol);end图像处理程序%测试霍夫变换clcclearcloseall%读取图像I =imread('circuit.tif');rotI=imrotate(I,80,'crop');%旋转33度，保持原图片大小fig1=imshow(rotI);%提取边BW=edge(rotI,'canny');figure,imshow(BW);%霍夫变换[H,theta,rho]=hough(BW);%计算二值图像的标准霍夫变换，H为霍夫变换矩阵，theta,rho为计算霍夫变换的角度和半径值figure,imshow(imadjust(mat2gray(H)),[],'XData',theta,'YData',rho,...'InitialMagnification','fit');xlabel('\theta(degrees)'),ylabel('\rho');axison,axisnormal,holdon;colormap(hot)%显示霍夫变换矩阵中的极值点P=houghpeaks(H,50,'threshold',ceil(0.3*max(H(:))));%从霍夫变换矩阵H中提取5个极值点x=theta(P(:,2));y=rho(P(:,1));plot(x