找到最佳价值

发布者Loren Shure.那8月3日，2021年

13.3.views (last 30 days) |0.Likes|13评论

            您是否曾经需要解决有当地最小值的优化问题？你用什么策略来解决它，试图找到“最好的”答案？今天我要谈论一个简单的策略，很容易提供
            全局优化工具箱
            。
           

解决一个简单的问题

            Or at least let's try. I have some data and I want to fit a particular form of a curve to it. First let's look at the pharmacokinetic data. Here's the reference:
            通过全局优化的非线性代数模型参数估计。
            电脑与化学工程，第22卷，1998年3月15日，第13页S213-S220 William R. Esposito，Christodoulos A. Floudas。
           

            数据是时间与浓度
           

              T.=[ 3.92, 7.93, 11.89, 23.90, 47.87, 71.91, 93.85, 117.84 ]
             

                  T.=
                  1×8
                 

                  3.9200 7.9300 11.8900 23.9000 47.8700 71.9100 93.8500 117.8400

              c = (0.163, 0.679, 0.679.那0.。3.8.8.那0.。18.3.那0.。125, 0.086, 0.0624 ]
             

                  C=
                  1×8
                 

                  0.1630 0.6790 0.6790 0.3880 0.1830 0.1250 0.0860 0.0624

            我喜欢看到数据，部分是为了确保我没有进入错误，部分是为了让整体系统的感受。事实上，让我们来看看数据。
           

              绘图（T，C，'o'的）
             

              xlabel('时间'的）
             

              ylabel('Concentration'的）
             

3个隔间模型

            As in the reference, we fit a 3 compartment model, sum of 3 decaying exponentials.
           

            c=b1 e(-b4 t)+b2 e(-b5 t)+b3 e(-b6 t)" style="vertical-align:-6px">
             
                C
               
                =
               
                  B.
                
                  1
                
                  E.
                
                   （
                  
                    -
                   
                      B.
                    
                      4.
                    
                    T.
                  
                   的）
                 
                +
               
                  B.
                
                  2
                
                  E.
                
                   （
                  
                      -
                     
                      B.
                    
                      5.
                    
                    T.
                  
                   的）
                 
                +
               
                  B.
                
                  3.
                
                  E.
                
                   （
                  
                      -
                     
                      B.
                    
                      6.
                    
                    T.
                  
                   的）

            我们可以表达那个模型
            匿名功能
            的
            T.
            （时间）和模型参数[
            B（1）B（2）... B（6）]
            。
           

              型号= @（b，t）b（1）* exp（-b（4）* t）+ b（2）* exp（-b（5）* t）+ b（3）* exp（-b（6）* t）
             

                 模型=function_handle具有值：@（b，t）b（1）* exp（-b（4）* t）+ b（2）* exp（-b（5）* t）+ b（3）* exp（-b（6）* t）

定义优化问题

            我们接下来定义要使用的优化问题
            problem-based formulation
            。这允许我们选择我们想要的解决者，提供数据和自然表达的约束和选项。
           

              问题= createOptimproblem（'lsqcurvefit'那......
             

              '客观的'那model,......
             

              'xdata'，t，'ydata'， C，......
             

              'x0'，一个（1,6），......
             

              'lb'，[-10 -10 -10 0 0]，......
             

              'UB'，[10 10 10 0.5 0.5 0.5]，......
             

              'options'，优化选择（'lsqcurvefit'那......
             

              'outputfcn'，@curvefittingplottates，......
             

              'Display'那'没有任何'））
             

                 problem =结构与字段：
                
                 目的：@（b，t）b（1）* exp（-b（4）* t）+ b（2）* exp（-b（5）* t）+ b（3）* exp（-b（6）* t）x0: [1 1 1 1 1 1] xdata: [3.9200 7.9300 11.8900 23.9000 47.8700 71.9100 93.8500 117.8400] ydata: [0.1630 0.6790 0.6790 0.3880 0.1830 0.1250 0.0860 0.0624] lb: [-10 -10 -10 0 0 0] ub: [10 10 10 0.5000 0.5000 0.5000] solver: 'lsqcurvefit' options: [1×1 optim.options.Lsqcurvefit]

解决这个问题

            First solve the problem directly once.
           

              B.=lsqcurvefit(problem)
             

                
                  B.=
                  1×6.
                 
                  0.1842 0.1836 0.1841 0.0172 0.0171 0.0171

            You'll notice that the model does not do a stellar job fitting the data or even following the shape of the data.
           

多层

            让我们看看我们是否可以通过从一堆不同的点开始做得更好。
           

              ms = multiStart;
             

              MS.Display ='iter';
             

              rng.默认
             

              数字
             

t

              [〜，Fval，ExitFlag，输出，解决方案万博 尤文图斯] =运行（MS，问题，50）
             

                运行本地本地本地一阶索引ExitFlag F（x）＃iter F计数最优性1 3 0.222 8 63 0.0006396 2 3 0.000154 21154 0.001864 3 3 0.009442 40 315 0.01989 4 3 1.462E-05 34 245 0.002586 5 3 1.454E-05 19 140 1.079E-05 6 3 1.475E-05 24 175 0.006883 7 3 0.009445 50 32 231 0.006853 1.495C-05 32 231 0.006853 1. 4.466C-05 32 231 0.00478 10 3 0.00944 80 567 0.01042 11 31.471E-05 40 287 0.005472 12 3 1.566E-05 24 175 0.001576 13 0.009439 24 3 0.0005121 14 3 0.029451 41 0.02935 15 3 0.009493 26 189 0.004837 3 1.476E-05 40 287 0.006352 17 3 1.494E-05 40287 0.008288 18 3 0.009446 62 441 0.01296 19 3 1.457E-05 22 161 0.001755 20 3 1.488E-05 53 378 0.007608 21 3 0.00944 32 22 0.006878 22 3 1.464E-05 24 175 0.003709 23 3 0.00949 43 308 0.02515 24 3 0.00944747 336 0.007942 25 3 1.455E-05 23 168 0.001621 22 0.009442 32 231 0.01328 23 1.479E-05 40 287 0.004 3 1.479E-05 4.4790-05 18 13. 0.09878 23 0. 0.0009721 30 3 1. 4.455E-05 42 301 0.001122 31 3 0.009441 47336 0.01537 32 3 0.009451 47 336 0.008942 33 3 0.0001729 14 105 0.0003276 34 3 0.009442 40 0.01062 32 31 0.01062 35 3 3. 0.09442 4 0.0109422 4 0.01452-05 26 189 0.009896 32 0.009896 32 0.009896 32 0.009458 39 280 0.02208 391824 175 7.815E-08 39 3 0.009441 60 427 0.0107 40 3 1.472E-05 34 245 0.002981 41 3 1.503E-05 22 161 0.00585 42 3 0.00952 15 112 0.008492 32 0.009439 21 0.009439 21 154 0.000769 42 3.462E-05 64 4550.0005576 45 3 0.009439 176 0.0001567 46 3 1.471E-05 30 217 0.0011973 47 3.009444 323 0.02022 48 3 1.474E-05 24 175 0.004799 49 3 1.522E-05 42 301 0.008228

                
               

                50 3 0.009445 40 287 0.02166 MultiStart完成了所有起点的运行。所有50个本地求解器运行融合，并使用正本地求解器退出标志。

               FVAL = 1.4540E-05

               EXITFLAG = 1

                 输出=结构与字段：
                

                 funccount：12726 localsolvertotal：50 localsolversuccess：50 localsolverincomplete：0 localsolvernosolution：0消息：'multiStart从所有起始点完成了运行.LALL 50本地解算器运行融合了正本地求解器退出标志。'

                 万博 尤文图斯解决方案=
                 1×50对象
                

                 
                  
                    
                    1
                    2
                    3.
                    4.
                    5.
                    6.
                    7.
                    8.
                    9.
                    10.
                    11.
                    12.
                    13.
                    14.
                    15.
                    16.
                    17.
                    18.
                    19.
                    20.
                    21.
                    22.
                    23.
                    24.
                    25.
                    26.
                    27.
                    28.
                    29.
                    30.
                   

                  
                    1
                    1×1 globaloptimsolution
                    1×1 globaloptimsolution
                    1×1 globaloptimsolution
                    1×1 globaloptimsolution
                    1×1 globaloptimsolution
                    1×1 globaloptimsolution
                    1×1 globaloptimsolution
                    1×1 globaloptimsolution
                    1×1 globaloptimsolution
                    1×1 globaloptimsolution
                    1×1 globaloptimsolution
                    1×1 globaloptimsolution
                    1×1 globaloptimsolution
                    1×1 globaloptimsolution
                    1×1 globaloptimsolution
                    1×1 globaloptimsolution
                    1×1 globaloptimsolution
                    1×1 globaloptimsolution
                    1×1 globaloptimsolution
                    1×1 globaloptimsolution
                    1×1 globaloptimsolution
                    1×1 globaloptimsolution
                    1×1 globaloptimsolution
                    1×1 globaloptimsolution
                    1×1 globaloptimsolution
                    1×1 globaloptimsolution
                    1×1 globaloptimsolution
                    1×1 globaloptimsolution
                    1×1 globaloptimsolution
                    1×1 globaloptimsolution
                   

                 
                

              serialtime = toc;
             

可视化最佳解决方案

            第50个解决方案，这是上面绘制的，不一定是最好的。幸运的是，
            多层
            从最坏的情况下订购解决万博 尤文图斯方案。所以我们只需要看第一个。
           

              CurveFittingPlotiterates（万博 尤文图斯解决方案）
             

            您现在可以看到第50架不是最好的解决方案，因为这一决赛的平均方形错误显示出更好的10倍。
           

多层with Parallel Computing

            我现在可以看看我们是否可以在本地使用所有4个核心作为平行工人提高性能。
           

              MS.USEPLASTER = TRUE;
             

              gcp;
             

              Tic;
             

              rng.默认
             

              运行（MS，问题，50）;
             

                Running the local solvers in parallel. Run Local Local Local Local First-order Index exitflag f(x) # iter F-count optimality 1 3 0.222 8 63 0.0006396 10 3 0.00944 80 567 0.01042 9 3 1.466e-05 35 252 0.00478 2 3 0.000154 21 154 0.001864 16 3 1.476e-05 40 287 0.006352 15 3 0.009493 26 189 0.004837 14 3 0.009451 41 294 0.02935 3 3 0.009442 44 315 0.01989 22 3 1.464e-05 24 175 0.003709 21 3 0.00944 37 266 0.006878 20 3 1.488e-05 53 378 0.007608 4 3 1.462e-05 34 245 0.002586 8 3 1.495e-05 32 231 0.006853 7 3 0.009445 50 357 0.002266 6 3 1.475e-05 24 175 0.006883 5 3 1.454e-05 19 140 1.079e-05 28 3 1.479e-05 18 133 0.006878 27 3 1.479e-05 40 287 0.004821 26 3 0.009442 32 231 0.01328 25 3 1.455e-05 23 168 0.001621 13 3 0.009439 24 175 0.0005121 12 3 1.566e-05 24 175 0.001576 11 3 1.471e-05 40 287 0.005472 34 3 0.009442 44 315 0.01062 33 3 0.0001729 14 105 0.0003276 32 3 0.009451 47 336 0.008942 19 3 1.457e-05 22 161 0.001755 18 3 0.009446 62 441 0.01296 17 3 1.494e-05 40 287 0.008288 24 3 0.009447 47 336 0.007942 23 3 0.009449 43 308 0.02515 31 3 0.009441 47 336 0.01537 30 3 1.455e-05 42 301 0.001122 29 3 1.456e-05 72 511 0.0009721 42 3 0.00952 15 112 0.008492 45 3 0.009439 17 126 0.0001567 37 3 0.009458 39 280 0.02208 36 3 1.509e-05 26 189 0.009896 35 3 0.0001751 21 154 7.71e-05 40 3 1.472e-05 34 245 0.002981 39 3 0.009441 60 427 0.0107 38 1 1.454e-05 24 175 7.815e-08 41 3 1.503e-05 22 161 0.00585 49 3 1.522e-05 42 301 0.008228 47 3 0.009444 38 273 0.02022 44 3 1.462e-05 64 455 0.0005576 43 3 0.009439 21 154 0.000769 48 3 1.474e-05 24 175 0.004799 50 3 0.009445 40 287 0.02166 46 3 1.471e-05 30 217 0.001973 MultiStart completed the runs from all start points. All 50 local solver runs converged with a positive local solver exit flag.

              parallelTime = toc;
             

Calculate Speedup

            由于池启动时间，速度可能不显而易见。自从我早先开始，我可以看到体面的加速。
           

              speedup = serialtime / parartateTime
             

               Speedup = 2.5014

Do You Have Problems Where the Solution is Sensitive to the Starting Point?

            Tell us
            关于您对解决方案空间的探索到您的问题。如果解决方案对您开始的位置敏感，则可能会考虑使用
            多层
            and other techniques from the Global Optimization Toolbox.
           

            版权所有2021 MathWorks，Inc。
           

附录

            这是绘制迭代的代码。
           

              dbtype.曲线特征普罗特
             

                1 = curvefittingPlotIterates (x, optim功能停止Values,state) 2 % Output function that plots the iterates of the optimization algorithm. 3 4 % Copyright 2010 The MathWorks, Inc. 5 6 persistent x0 r; 7 if nargin == 1 8 showPlot(x(1).X,x(1).X0{:},x(1).Fval) 9 else 10 switch state 11 case 'init' % store initial point for later use 12 x0 = x; 13 case 'done' 14 if ~(optimValues.iteration == 0) 15 % After optimization, display solution in plot title 16 r = optimValues.resnorm; 17 showPlot(x,x0,r) 18 end 19 end 20 end 21 if nargout > 0 22 stop = false; 23 clear function 24 end 25 end 26 27 function showPlot(b,b0,r) 28 f = @(b,x) b(1)*exp(-b(4).*x) + b(2).*exp(-b(5).*x) +... 29 b(3).*exp(-b(6).*x); 30 31 persistent h ha 32 if isempty(h) || ~isvalid(h) 33 x = [ 3.92, 7.93, 11.89, 23.90, 47.87, 71.91, 93.85, 117.84 ]; 34 y = [ 0.163, 0.679, 0.679, 0.388, 0.183, 0.125, 0.086, 0.0624 ]; 35 plot(x,y,'o'); 36 xlabel('t') 37 ylabel('c') 38 title('c=b_1e^{-b_4t}+b_2e^{-b_5t}+b_3e^{-b_6t}') 39 axis([0 120 0 0.8]); 40 h = line(3:120,f(b,3:120),'Color','r','Tag','PlotIterates'); 41 42 else 43 set(h,'YData',f(b,get(h,'XData'))); 44 end 45 s = sprintf('Starting Value Fitted Value\n\n'); 46 47 for i = 1:length(b) 48 s = [s, sprintf('b(%d): % 2.4f b(%d): % 2.4f\n',i,b0(i),i,b(i))]; 49 end 50 s = [s,sprintf('\nMSE = %2.4e',r)]; 51 52 if isempty(ha) || ~isvalid(ha) 53 % Create textbox 54 ha = annotation(gcf,'textbox',... 55 [0.5 0.5 0.31 0.32],... 56 'String',s,... 57 'FitBoxToText','on',... 58 'Tag','CoeffDisplay'); 59 end 60 ha.String = s; 61 drawnow 62 63 end