使用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*pi*2.5;Amp = 2.2;相位= -[0;0.54;2.07);numSources = 3;身高= 3;%所有点源对齐于[x_i,y_i,z]xcods = [2.4112 0.2064 1.6787];ycods = [0.3957 0.3927 0.9877];zcods = height*ones(numSources,1);起源= [xcods ycods zcods];

想象干涉模式

现在让我们在z = 0平面上看到干涉图的一个切片。

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

通过匿名函数传递附加参数:waveIntensity_x = @(x).．.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})²+ (y - y_{c1})²\le r_1²$

$(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)，域，150);持有在%绘图约束G1 = @(x,y) (x- xcoordinates(1))。^2 + (y- ycoordinates (1))^2 - 9;G2 = @(x,y) (x- xcoordinates(2))。^2 + (y- ycoordinates (2))^2 - 9;G3 = @(x,y) (x- xcoordinates(3))。^2 + (y- ycoordinates (3))^2 - 9;H1 = ezplot(g1，域);h1。Color = [0.8 0.7 0.1];%的黄色h1。线宽= 1.5;H2 = ezplot(g2，域);h2。Color = [0.3 0.7 0.5];%绿色h2。线宽= 1.5;H3 = ezplot(g3，域);h3。Color = [0.4 0.4 0.6];%的蓝色h3。线宽= 1.5;%绘图边界Lb = [-0.5 -2];Ub = [3.5 3];行([lb(1) lb(1)]，[lb(2) ub(2)]，“线型”，“——”) line([ub(1) ub(1)]，[lb(2) ub(2)]，“线型”，“——”) line([lb(1) ub(1)]，[lb(2) lb(2)]，“线型”，“——”) line([lb(1) ub(1)]，[ub(2) ub(2)]，“线型”，“——”)标题(“带约束边界的模式轮廓”）

用本地求解器设置和解决问题

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

让我们建立一个问题结构来描述我们的优化问题。我们想要最大化强度函数，所以我们对返回的值求反waveIntensity．让我们任意选择一个在可行域附近的起点。

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

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

找到满足约束条件的局部最小值。优化完成是因为目标函数在可行方向上不递减，在最优性容差值范围内，约束条件满足在约束容差值范围内。Xlocal = -0.5000 fvallocal = -1.4438

现在，让我们看看我们是如何显示的结果fmincon在等高线图中。请注意,fmincon没有达到全局最大值，这也是在图上标注的。注意，我们只画出在解处活跃的边界。

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

Xmax = 1.2500 0.4450

使用`GlobalSearch`而且`MultiStart`

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

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

构造一个GlobalSearch对象gs = GlobalSearch;根据GlobalSearch属性构造一个MultiStart对象ms = MultiStart;rng (4“旋风”）再现率%%运行GlobalSearch抽搐;[xgs，~，~，~，solsgs] = run(gs,problem);toc xg使用15个随机生成的点运行MultiStart抽搐;[xms，~，~，~，solsms] = run(ms,problem,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结果xGS = 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我们能够在这次运行中定位全局最大值。

当对目标函数或约束条件知之甚少时，在实践中很难找到严密的边界。不过，一般来说，我们可以猜测出一个合理的区域，我们希望在其中限制起始点集。出于说明目的，让我们放宽边界，定义一个更大的区域，在其中生成起点并重新尝试求解器。

放宽边界以分散起始点。问题。Lb = -5*ones(2,1);问题。Ub = 5*ones(2,1);%运行GlobalSearch抽搐;[xgs，~，~，~，solsgs] = run(gs,problem);toc xg使用15个随机生成的点运行MultiStart抽搐;[xms，~，~，~，solsms] = run(ms,problem,15);toc xms

GlobalSearch停止了，因为它分析了所有的试验点。所有4个本地求解器运行都收敛于一个积极的本地求解器退出标志。运行时间为0.173760秒。MultiStart完成了所有起点的运行。所有15个本地求解器运行都收敛于一个积极的本地求解器退出标志。运行时间为0.134150秒。XMS = 2.4947 -0.1439

显示轮廓身材轮廓(X,Y,Z) hold住在创建显示xglobal位置的文本箭头注释(“textarrow”，[0.44 0.50]，[0.63 0.58]，“TextEdgeColor”，[0 0 0]，.．.“TextBackgroundColor”，[11 11 1]，“字形大小”12“字符串”, {“全球最大”}）;轴([-1 3.75 -3 3])使用“*”标记绘制GlobalSearch结果xGS = 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”，“最好的女士”，“位置”，“最佳”)标题(“GlobalSearch和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抽搐;XGS = run(gs,problem);toc xg

GlobalSearch停止了，因为它分析了所有的试验点。所有10个本地求解器运行都收敛于一个积极的本地求解器退出标志。运行时间为0.242955秒。XGS = 1.2592 0.4284

利用`MultiStart`并行能力

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

设置MultiStart的UseParallel属性ms = MultiStart(ms，“UseParallel”,真正的);试一试demoOpenedPool = false;%如果不存在，请创建一个并行池%(需要并行计算工具箱)如果Max (size(gcp)) == 0%，如果没有池parpool demoOpenedPool = true;结束抓我警告消息(“globaloptim: globaloptimdemos: opticalInterferenceDemo: noPCT”));结束运行求解器抽搐;XMS = run(ms,problem,100);toc xms如果demoOpenedPool如果在本例中创建了池，请确保删除池删除(gcp)删除存储池结束

MultiStart从所有起始点完成运行。所有100个本地求解器运行都收敛于一个正的本地求解器退出标志。运行时间为0.956671秒。XMS = 1.2592 0.4284

客观约束和非线性约束

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

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

函数[c,ceq] = apertureConstraint(x, xcoordinates, ycoordinates)光干涉演示的孔径约束函数。版权所有The MathWorks, Inc.Ceq = [];C = (x(1) - xcoordinates)。^2 + (x(2) - ycoordinates)^2 - 9;

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