罗兰谈MATLAB的艺术

将想法转化为MATLAB

寻找最优值

你曾经需要解决一个有局部极小值的优化问题吗?你用什么策略来解决它,试图找到“最佳”答案?今天我要讲一个简单的策略,在 全局优化工具箱

解决简单问题

或者至少让我们试一试。我有一些数据,我想用曲线的特定形式来拟合它。首先让我们看一下药代动力学数据。参考资料如下: 基于全局优化的非线性代数模型参数估计。 计算机与化学工程,第22卷,增刊1,1998年3月15日,第S213-S220页
数据是时间和浓度
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.388, 0.183, 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”
包含(“时间”
ylabel (“浓度”

3车厢模型

正如在参考文献中,我们拟合一个3室模型,3个衰减指数的和。
c=b1e(-b4t)+b2e(-b5t)+b3e(-b6t)" 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(2)…b (6)]
模型= @ (b, t) b (1) * exp (t - b (4) *) + b (2) * exp (t - b (5) *) + b (3) * exp (t - b (6) *)
模型=Function_handle with value:@ (b, b t) (1) * exp (t - b (4) *) + b (2) * exp (t - b (5) *) + b (3) * exp (t - b (6) *)

定义优化问题

方法定义要解决的优化问题 具体问题具体分析公式 .这允许我们选择我们想要的求解器,提供数据,并自然地表达约束和选项。
问题= createOptimProblem(“lsqcurvefit”...
“目标”、模型...
“xdata”t“ydata”c...
“x0”的(6),...
“磅”, [-10 -10 -10 0 0 0],...
乌兰巴托的, [10 10 10 0.5 0.5 0.5],...
“选项”optimoptions (“lsqcurvefit”...
“OutputFcn”@curvefittingPlotIterates,...
“显示”“没有”))
问题=带字段的结构:
目标:@(b,t)b(1)*exp(-b(4)*t)+b(2)*exp(-b(5)*t)+b(3)*exp(-b(6)*t) x0: [1 11 11 11 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 000 0] ub: [10 10 10 0.5000 0.5000 0.5000 0.5000] solver: 'lsqcurvefit' options: [1×1 optim.options.Lsqcurvefit]

解决问题

首先直接解决问题一次。
B = lsqcurvefit(问题)
b = 1×6
0.1842 0.1836 0.1841 0.0172 0.0171 0.0171
您将注意到,该模型在拟合数据或遵循数据形状方面做得并不出色。

MultiStart

我们看看从不同的点开始是否能做得更好。
ms = MultiStart;
ms.Display =“通路”
rng默认的
数字
抽搐
[~,fval,exitflag,output,金宝搏官方网站solutions] = run(ms, problem, 50)
运行当地当地当地一阶指数exitflag f (x) # iter F-count最优1 3 0.222 8 0.000154 21 154 0.001864 63 0.0006396 - 2 3 3 3 0.009442 44 315 0.01989 4 3 e-05 34 245 0.002586 1.462 1.454 e-05 19 140 1.079 e-05 6 3 1.475 e-05 24 175 0.006883 7 3 0.009445 357 0.002266 8 3 1.495 e-05 32 231 1.466 0.006853 9 3 e-05 35 252 0.00478 567 0.00944 80 0.01042 11 3 1.471 e-05 175 0.001576 287 1.566 0.005472 12 3 e-05 24 13 3 0.009439 24 175 0.0005121 14 3 0.009451 41 294 0.02935 153 0.009493 26 189 0.004837 16 3 1.476 e-05 287 0.008288 287 1.494 0.006352 17 3 e-05 40 18 62 441 0.009446 1.457 0.01296 19 3 e-05 22 161 0.001755 378 0.007608 1.488 e-05 53 21 22 3 0.00944 37 266 0.006878 175 0.003709 1.464 e-05 24 23 3 0.009449 43 308 0.02515 24 336 0.007942 0.009447 47 25 3 1.455 e-05 23 168 0.009442 0.001621 26 3 32 231 0.01328 27日3 1.479 e-05 40 287 0.004821 28日3 1.479 e-05 18 133 0.006878 1.456 e-05 72 29 511 0.0009721 301 0.001122 1.455 e-05 42 31 3 0.009441 473.3.60.0153.7323.0.009451 47 336 0.008942 33 3 0.0001729 14 105 0.0003276 34 3 0.009442 44 315 0.01062 35 3 0.0001751 21 154 7.71e-05 36 3 1.509e-05 26 189 0.009896 37 3 0.009458 39 280 0.02208 38 1 1.454e-05 24 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 43 3 0.009439 21 154 0.000769 44 3 1.462e-05 64 455 0.0005576 45 3 0.009439 17 126 0.0001567 46 3 1.471e-05 30 217 0.001973 47 3 0.009444 38 273 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 message: 'MultiStart从所有起始点完成运行。(()所有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个解,也就是上面画的,不一定是最好的解。幸运的是, MultiStart 从最佳到最差排序。金宝搏官方网站我们只需要看第一个。
curvefittingPlotIterates(金宝搏官方网站解决方案)
现在你可以看到,第50个不是最好的解决方案,因为最后一个显示的均方误差超过10倍。

MultiStart与并行计算

现在,我将看看是否可以将我的所有4个核心作为本地并行工作来提高性能。
ms.UseParallel = true;
gcp;
抽搐;
rng默认的
运行(ms, problem, 50);
并行运行局部求解器。运行当地当地当地一阶指数exitflag f (x) # iter F-count最优1 3 8 0.222 63 0.0006396 567 0.00944 80 0.01042 252 0.00478 1.466 e-05 35 2 3 154 21 0.000154 1.476 0.001864 16 3 e-05 189 0.004837 287 0.009493 0.006352 15 3 26 14 3 0.009451 41 294 0.02935 315 0.01989 0.009442 44 22 3 1.464 e-05 24 175 0.003709 21 3 0.00944 266 0.006878 378 0.007608 1.488 e-05 53 4 3 1.462 e-05 34 245 0.002586 8 3 1.495 e-05 32 231 0.009445 357 0.002266 0.006853 7 3 6 31.475 e-05 24 175 0.006883 140 1.079 1.454 e-05 19 e-05 28 133 1.479 0.006878 27 3 1.479 e-05 18 e-05 231 0.01328 287 0.009442 0.004821 26 3 32 25 3 1.455 e-05 23 24 175 168 0.001621 13 3 0.009439 1.566 0.0005121 12 3 e-05 24 175 1.471 0.001576 11 3 e-05 40 287 0.009442 0.005472 34 3 44 315 0.01062 33 3 0.0001729 336 0.008942 105 0.009451 0.0003276 32 3 47 19 3 18 1.457 e-05 22 161 0.001755 0.009446 0.008288 62 441 1.494 0.01296 17 3 e-05 40 287 24 336 0.007942 0.009447 47 23 3 0.009449 4330.80.025153.13. 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;

计算加速

由于池启动时间的关系,在第二次运行之前,速度可能不会明显提高。因为我开始得更早,我看到了不错的速度。
speedup = serialTime/parallelTime
加速= 2.5014

你是否遇到了解决方案对起点敏感的问题?

告诉我们 关于你对问题解决空间的探索。如果解决方案对开始的位置很敏感,则可以考虑使用 MultiStart 以及全局优化工具箱中的其他技术。
The MathWorks, Inc.版权所有

附录

下面是绘制迭代的代码。
dbtypecurvefittingPlotIterates
1函数stop = curvefittingPlotIterates(x,optimValues,state) 2%绘制优化算法迭代次数的输出函数。3 4%版权所有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开关状态11 case 'init' %存储初始点以供以后使用12 x0 = x;13 case 'done' 14 if ~(迭代== 0)15%优化后,在图标题16中显示解r = optimValues.resnorm;17 showPlot(x,x0,r) 18 end 19 end 20 end 21 if nargout > 0 22 stop = false;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
|
  • 打印
  • 发送电子邮件

评论

如欲留言,请点击在这里登录您的MathWorks帐户或创建一个新帐户。