利用GlobalSearch和MultiStart最大化单色偏振光干涉图样

打开脚本

这个例子展示了如何使用这些函数GlobalSearch而且MultiStart．

简介

这个例子特别说明了全局优化工具箱是如何工作的GlobalSearch而且MultiStart，可以帮助定位电磁干扰图的最大值。为了简化建模，图案产生于从点源扩散出来的单色偏振光。

在点x和时间t处，由源i在偏振方向上测量的电场为

$$ E_i = \压裂{ai} {d_i (x)} \罪(\ phi_i + \ω(t - d1 (x) / c)), $$

在哪里 $\ phi_i美元$ 源的相位在时间为0时吗我美元，美元加元是光速， $ω\美元$ 是光的频率， ai美元是源的振幅吗, d_i (x)美元是到源的距离吗来 x美元．

对于一个不动点 x美元光的强度是净电场平方的时间平均值。净电场是所有源产生的电场之和。平均时间只取决于电场的大小和相对相位 x美元．为了计算净电场，用相量法把各个贡献相加。对于相量，每个源贡献一个矢量。矢量的长度是振幅除以到源的距离，矢量的角度， $\phi_i - \omega d_i(x)/c$ 是该点的相位。

对于本例，我们定义了三个具有相同频率的点源( $ω\美元$ )和振幅( 一个美元 )，但初始阶段不同( $\ phi_i美元$ )．我们把这些源排列在一个固定的平面上。

%频率与峰值的数量成正比relFreqConst = 2 *π* 2.5;amp = 2.2;阶段= - [0;0.54;2.07);numSources = 3;身高= 3;%所有点源对齐在[x_i,y_i,z]Xcoords = [2.4112 0.2064 1.6787];Ycoords = [0.3957 0.3927 0.9877];zcoords =身高*的(numSources, 1);Origins = [xcoords ycoords zcoords];

干涉图样形象化

现在让我们想象一下z = 0平面上干涉图样的切片。

从下图中可以看出，有许多波峰和波谷表示建设性干扰和破坏性干扰。

通过匿名函数传递附加参数:waveIntensity_x = @(x) waveIntensity(x,amp,phase，.．.relFreqConst numSources,起源);%生成网格(X, Y) = meshgrid (4:0.035:4 4:0.035:4);计算网格上的强度Z = arrayfun(@(x,y) waveIntensity_x([x,y])， x,y);画出曲面和等高线图冲浪(X, Y, Z,“EdgeColor”，“没有”)包含(“x”) ylabel (“y”) zlabel (“强度”）

提出优化问题

我们感兴趣的是这个波的强度达到最高点的位置。

波浪强度( 我美元 )在远离源时下降，与 1美元/ d1 (x)美元．因此，让我们通过向问题添加约束来限制可行解决方案的空间。

如果我们用一个光圈来限制光源的曝光，那么我们就可以预期光圈在观测平面上投影的交点处会得到最大值。我们通过将搜索限制在以每个光源为中心的圆形区域来模拟孔径的效应。

我们还通过给问题增加边界来限制解空间。尽管这些边界可能是多余的(考虑到非线性约束)，但它们是有用的，因为它们限制了起始点生成的范围(参见GlobalSearch和MultiStart是如何工作的的更多信息)。

现在我们的问题变成了:

$\max_{x,y} I(x,y)$

受

$(x - x_{c1})^2 + (y - y_{c1})^2 \le r_1^2$

$(x - x_{c2})^2 + (y - y_{c2})^2 \le r_2^2$

$(x - x_{c3})^2 + (y - y_{c3})^2 \le r_3^2$

$-0.5 \leq x \leq 3.5$

$-2 \leq y \leq 3$

在哪里 $y_(美元间{cn}, {cn})美元$ 而且 r_n美元是坐标和孔径半径吗 $n ^ {th}$ 分别点光源。每个光源的光圈半径为3。给定的边界包含可行区域。

我们的目标( 我(x, y)美元 )和非线性约束函数在单独的MATLAB®文件中定义，waveIntensity.m而且apertureConstraint.m，它们分别列在本例的末尾。

可视化与约束

现在让我们用非线性约束边界的叠加来可视化干涉图的轮廓。可行区域是三个圆(黄、绿、蓝)相交的内部。变量的边界由虚线框表示。

形象化干涉曲面的轮廓Domain = [-3 5.5 -4 5];图;ezcontour(@(X,Y) arrayfun(@(X,Y) waveIntensity_x([X,Y])，X,Y)，domain,150);持有在%的阴谋约束G1 = @(x,y) (x- xcods(1))。^ 2 + (y-ycoords(1))。^ 2 - 9;G2 = @(x,y) (x- xcoors(2))。^ 2 + (y-ycoords(2))。^ 2 - 9;G3 = @(x,y) (x- xcoors(3))。^ 2 + (y-ycoords(3))。^ 2 - 9;h1 = ezplot (g1、领域);h1。颜色= [0.8 0.7 0.1];%的黄色h1。线宽= 1.5;h2 = ezplot (g2、领域);h2。颜色= [0.3 0.7 0.5];%绿色h2。线宽= 1.5;h3 = ezplot (g3、领域);h3。颜色= [0.4 0.4 0.6];%的蓝色h3。线宽= 1.5;%画界限Lb = [-0.5 -2];Ub = [3.5 3];线([磅(1)磅(1)]、[磅(2)乌兰巴托(2)),“线型”，“——”)线([ub(1) ub(1)]，[lb(2) ub(2)]，“线型”，“——”)线([lb(1) ub(1)]，[lb(2) lb(2)]，“线型”，“——”)线([lb(1) ub(1)]，[ub(2) ub(2)]，“线型”，“——”)标题(“带约束边界的图案轮廓”）

用本地求解器建立和解决问题

给定非线性约束，我们需要一个有约束的非线性求解器，即fmincon．

让我们建立一个问题结构来描述我们的优化问题。我们想要使强度函数最大化，所以我们对返回的值进行否定waveIntensity．我们随便选一个起点刚好在可行区域附近。

对于这个小问题，我们将使用fmincon的SQP算法。

通过匿名函数传递附加参数:apertureConstraint_x = @(x) apertureConstraint(x,xcoords,ycoords);设置fmincon的选项X0 = [3 -1];选择= optimoptions (“fmincon”，“算法”，“sqp”）;问题= createOptimProblem (“fmincon”，“目标”，.．.@ (x) -waveIntensity_x (x)“x0”x0,“磅”磅,乌兰巴托的乌兰巴托,.．.“nonlcon”apertureConstraint_x,“选项”、选择);%叫fmincon[xlocal, fvallocal] = fmincon(问题)

找到了满足约束条件的局部极小值。由于目标函数在可行方向上不减少，优化完成，在最优性公差的值内，约束满足在约束公差的值内。Xlocal = -0.5000 0.4945 fvallocal = -1.4438

现在，让我们看看我们是怎么做的fmincon在等高线图中。请注意,fmincon没有达到全局最大值，这在图上也有注解。注意，我们只绘制解处活动的边界。

[~, maxIdx] = max (Z (:));xmax = [X(maxIdx)，Y(maxIdx)]图形轮廓(X,Y,Z)保持在%显示范围线([磅(1)磅(1)]、[磅(2)乌兰巴托(2)),“线型”，“——”）创建显示xlocal位置的文本箭头注释(“textarrow”(0.25 - 0.21), (0.86 - 0.60),“TextEdgeColor”(0 0 0),.．.“TextBackgroundColor”(1 1 1),“字形大小”11“字符串”, {“单一运行结果”});创建显示xglobal位置的文本箭头注释(“textarrow”(0.44 - 0.50), (0.63 - 0.58),“TextEdgeColor”(0 0 0),.．.“TextBackgroundColor”(1 1 1),“字形大小”12“字符串”, {“全球最大”});轴([-1 3.75 -3 3])

Xmax = 1.2500 0.4450

使用`GlobalSearch`而且`MultiStart`

给定一个任意的初始猜测，fmincon被卡在了附近的局部最大值。特别是全局优化工具箱求解器GlobalSearch而且MultiStart，让我们有更好的机会找到全局最大值，因为他们会尝试fmincon从多个生成的初始点(或者我们自己的自定义点，如果我们选择的话)。

我们的问题已经在问题结构，现在我们构造求解器对象并运行它们。的第一个输出运行是找到最佳结果的位置。

构造一个GlobalSearch对象gs = GlobalSearch;基于我们的GlobalSearch属性构造一个MultiStart对象女士= MultiStart;rng (4“旋风”）%的再现性% GlobalSearch运行抽搐;[xg, ~, ~, ~, solsgs] =运行(gs、问题);toc xg运行MultiStart与15个随机生成的点抽搐;[xms, ~, ~, ~, solsms] = (ms,问题,15);toc xms

GlobalSearch停止了，因为它分析了所有的试验点。所有14次本地求解器运行都带有一个正的本地求解器出口标志。运行时间为0.229525秒。MultiStart从所有起始点完成运行。所有15次本地求解器运行都汇聚了一个正的本地求解器出口标志。运行时间为0.109984秒。XMS = 1.2592 0.4284

检查结果

让我们检查一下两个求解器返回的结果。需要注意的重要一点是，根据为每个求解器创建的随机起点，结果会有所不同。对这个示例的另一次运行可能会得到不同的结果。最佳结果的坐标xg而且xms打印到命令行。我们将显示返回的惟一结果GlobalSearch而且MultiStart并从接近全局解决方案的角度突出每个求解器的最佳结果。

每个求解器的第五个输出是一个包含不同最小值(在本例中是极大值)的向量。我们画出结果的(x,y)对，solsgs而且solsms与我们之前用过的等高线图相反。

使用“*”标记绘制GlobalSearch结果xg = cell2mat ({solsgs(:)。x}”);散射(xg (: 1) xg (:, 2),‘*’，“MarkerEdgeColor”(0 0 1),“线宽”, 1.25)使用圆形标记绘制MultiStart结果xMS = cell2mat ({solsms(:)。x}”);散射(xMS (: 1), xMS (:, 2),“o”，“MarkerEdgeColor”(0 0 0),“线宽”传说,1.25)(“强度”，“约束”，“GlobalSearch”，“MultiStart”，“位置”，“最佳”)标题(“GlobalSearch和MultiStart结果”）

放松的界限

问题的边界很紧，两者都有GlobalSearch而且MultiStart我们能够在这次运行中找到全局最大值。

当对目标函数或约束条件了解不多时，在实践中很难找到较紧的边界。一般来说，我们可以猜测出一个合理的区域，在这个区域内我们可以限制起始点的集合。出于演示的目的，让我们放宽边界来定义一个更大的区域，在这个区域中生成起始点并重新尝试求解器。

放宽边界以分散起始点。问题。磅= 5 * 1 (2,1);问题。乌兰巴托= 5 * 1 (2,1);% GlobalSearch运行抽搐;[xg, ~, ~, ~, solsgs] =运行(gs、问题);toc xg运行MultiStart与15个随机生成的点抽搐;[xms, ~, ~, ~, solsms] = (ms,问题,15);toc xms

GlobalSearch停止了，因为它分析了所有的试验点。所有4次本地解算器运行都汇聚了一个正的本地解算器出口标志。运行时间为0.173760秒。MultiStart从所有起始点完成运行。所有15次本地求解器运行都汇聚了一个正的本地求解器出口标志。运行时间为0.134150秒。XMS = 2.4947 -0.1439

显示轮廓图轮廓(X, Y, Z)在创建显示xglobal位置的文本箭头注释(“textarrow”(0.44 - 0.50), (0.63 - 0.58),“TextEdgeColor”(0 0 0),.．.“TextBackgroundColor”(1 1 1),“字形大小”12“字符串”, {“全球最大”});轴([-1 3.75 -3 3])使用“*”标记绘制GlobalSearch结果xg = cell2mat ({solsgs(:)。x}”);散射(xg (: 1) xg (:, 2),‘*’，“MarkerEdgeColor”(0 0 1),“线宽”, 1.25)使用圆形标记绘制MultiStart结果xMS = cell2mat ({solsms(:)。x}”);散射(xMS (: 1), xMS (:, 2),“o”，“MarkerEdgeColor”(0 0 0),“线宽”, 1.25)%突出显示每一项的最佳结果:% GlobalSearch结果为红色，MultiStart结果为蓝色情节(xg (1) xg (2),“某人”，“MarkerSize”12“MarkerFaceColor”，[1 0 0]) plot(xms(1)，xms(2)，“某人”，“MarkerSize”12“MarkerFaceColor”，[0 0 1])传说(“强度”，“GlobalSearch”，“MultiStart”，“最好的GS”，“最好的女士”，“位置”，“最佳”)标题(“全球搜索和MultiStart结果与放松的边界”）

最好的结果GlobalSearch所示的红色方块和最佳结果来自MultiStart由蓝色方块显示。

调优`GlobalSearch`参数

注意，在这次运行中，由于边界定义了更大的区域，两个求解器都无法确定最大强度点。我们可以用几种方法来克服这个问题。首先,我们检查GlobalSearch．

请注意,GlobalSearch只跑了fmincon几次。为了增加找到全局最大值的机会，我们希望运行更多的点。为了将起始点集限制为最可能找到全局最大值的候选点，我们将指示每个求解器忽略不满足约束的起始点StartPointsToRun财产bounds-ineqs．此外，我们将设置MaxWaitCycle而且BasinRadiusFactor属性,这样GlobalSearch将能够快速识别狭窄的山峰。减少MaxWaitCycle原因GlobalSearch减小吸力半径BasinRadiusFactor比默认设置更频繁。

增加候选点总数，但过滤掉不可行点。gs = GlobalSearch (gs,“StartPointsToRun”，“bounds-ineqs”，.．.“MaxWaitCycle”3,“BasinRadiusFactor”, 0.3);% GlobalSearch运行抽搐;xg =运行(gs、问题);toc xg

GlobalSearch停止了，因为它分析了所有的试验点。所有10次本地求解器运行都汇聚了一个正的本地求解器出口标志。运行时间为0.242955秒。XGS = 1.2592 0.4284

利用`MultiStart`的并行功能

提高我们找到全局最大值的机会的一种强力方法是尝试更多的起始点。同样，这可能并不适用于所有情况。在我们的例子中，到目前为止我们只尝试了一小部分，运行时间也不是很长。所以，尝试更多的起点是合理的。为了加快我们的计算速度MultiStart并行，如果并行计算工具箱™是可用的。

设置MultiStart的UseParallel属性= MultiStart女士(女士,“UseParallel”,真正的);试一试demoOpenedPool = false;如果不存在，则创建一个并行池。%(需要并行计算工具箱)如果max(大小(gcp)) = = 0%，如果没有池parpool demoOpenedPool = true;结束抓我警告消息(“globaloptim: globaloptimdemos: opticalInterferenceDemo: noPCT”));结束运行求解程序抽搐;xms =运行(女士,100年问题);toc xms如果demoOpenedPool如果在本例中创建了池，请确保删除池删除(gcp)%删除池结束

MultiStart从所有起始点完成运行。所有100次本地求解器运行都汇聚了一个正的本地求解器出口标志。运行时间为0.956671秒。XMS = 1.2592 0.4284

客观约束与非线性约束

下面我们列出了定义优化问题的函数:

函数p = waveIntensity (x, amp,阶段,relFreqConst numSources,起源)光干涉演示的强度函数。版权所有2009 The MathWorks, Inc.d = distanceFromSource (x, numSources起源);ampVec =[总和(amp。/d .* cos(phase - d*relFreqConst));总和(amp。/d .* sin(phase - d*relFreqConst))];%强度为||AmpVec||^2p = ampVec ' * ampVec;

函数测查[c] = apertureConstraint (x, xcoords ycoords)孔径约束函数用于opalinterferencedemo。版权所有2009 The MathWorks, Inc.测查= [];C = (x(1) - xcoords)。^2 + (x(2) - ycoords)^ 2 - 9;

函数d = distanceFromSource (v, numSources起源)opticalInterferenceDemo的距离函数。版权所有2009 The MathWorks, Inc.d = 0 (numSources, 1);为k = 1:numSources d(k) = norm(origins(k，:) - [v 0]);结束