{"id":4418,"date":"2021-08-03T08:10:46","date_gmt":"2021-08-03T12:10:46","guid":{"rendered":"https:\/\/blogs.mathworks.com\/loren\/?p=4418"},"modified":"2021-08-05T18:19:14","modified_gmt":"2021-08-05T22:19:14","slug":"finding-the-optimal-value","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/loren\/2021\/08\/03\/finding-the-optimal-value\/","title":{"rendered":"Finding the Optimal Value"},"content":{"rendered":"
Have you ever needed to solve an optimization problem where there were local minima? What strategy do you use to solve it, trying to find the \"best\" answer? Today I'm going to talk about a simple strategy, readily available in the <\/span>Global Optimization Toolbox<\/span><\/a>.<\/span><\/div>

Solve a Simple Problem<\/span><\/h2>
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: <\/span>Parameter estimation in nonlinear algebraic models via global optimization.<\/span> Computers & Chemical Engineering, Volume 22, Supplement 1, 15 March 1998, Pages S213-S220 William R. Esposito, Christodoulos A. Floudas.<\/span><\/div>
The data are time vs. concentration<\/span><\/div>
t = [ 3.92, 7.93, 11.89, 23.90, 47.87, 71.91, 93.85, 117.84 ]<\/span><\/span><\/div>
t = <\/span>1\u00d78<\/span><\/div><\/div>
3.9200 7.9300 11.8900 23.9000 47.8700 71.9100 93.8500 117.8400\r\n<\/div>
<\/div>
<\/div><\/div><\/div><\/div><\/div><\/div>
c = [0.163, 0.679, 0.679, 0.388, 0.183, 0.125, 0.086, 0.0624 ]<\/span><\/span><\/div>
c = <\/span>1\u00d78<\/span><\/div><\/div>
0.1630 0.6790 0.6790 0.3880 0.1830 0.1250 0.0860 0.0624\r\n<\/div>
<\/div>
<\/div><\/div><\/div><\/div><\/div><\/div><\/div>
I like to see the data, in part to be sure I have no entry mistakes, and in part to get a feel for the overall system. In fact, let's visualize the data.<\/span><\/div>
plot(t,c,<\/span>'o'<\/span>)<\/span><\/span><\/div><\/div>
xlabel(<\/span>'Time'<\/span>)<\/span><\/span><\/div><\/div>
ylabel(<\/span>'Concentration'<\/span>)<\/span><\/span><\/div>
<\/div><\/div><\/div><\/div><\/div>

3 Compartment Model<\/span><\/h2>
As in the reference, we fit a 3 compartment model, sum of 3 decaying exponentials.<\/span><\/div>
c<\/mi>=<\/mo>b<\/mi><\/mrow>1<\/mn><\/mrow><\/msub>\u2009<\/mtext>e<\/mi><\/mrow>(<\/mo>-<\/mo>b<\/mi><\/mrow>4<\/mn><\/mrow><\/msub>\u2009<\/mtext>t<\/mi><\/mrow>)<\/mo><\/mrow><\/mrow><\/msup>+<\/mo>b<\/mi><\/mrow>2<\/mn><\/mrow><\/msub>\u2009<\/mtext>e<\/mi><\/mrow>(<\/mo>-<\/mo>b<\/mi><\/mrow>5<\/mn><\/mrow><\/msub>\u2009<\/mtext>t<\/mi><\/mrow>)<\/mo><\/mrow><\/mrow><\/msup>+<\/mo>b<\/mi><\/mrow>3<\/mn><\/mrow><\/msub>\u2009<\/mtext>e<\/mi><\/mrow>(<\/mo>-<\/mo>b<\/mi><\/mrow>6<\/mn><\/mrow><\/msub>\u2009<\/mtext>t<\/mi><\/mrow>)<\/mo><\/mrow><\/mrow><\/msup><\/mrow><\/math><\/span><\/div>
and we can express that model as an <\/span>anonymous function<\/span><\/a> of<\/span> t<\/span> (time) and the model parameters [<\/span>b(1) b(2) ... b(6)]<\/span>.<\/span><\/div>
model = @(b,t) b(1)*exp(-b(4)*t) + b(2)*exp(-b(5)*t) + b(3)*exp(-b(6)*t)<\/span><\/span><\/div>
model = function_handle with value:<\/span><\/span>\r\n @(b,t)b(1)*exp(-b(4)*t)+b(2)*exp(-b(5)*t)+b(3)*exp(-b(6)*t)<\/div><\/div><\/div><\/div><\/div><\/div>

Define Optimization Problem<\/span><\/h2>
We next define the optimization problem to solve using the <\/span>problem-based formulation<\/span><\/a>. This allows us to choose the solver we want, supply the data, and naturally express constraints and options.<\/span><\/div>
problem = createOptimProblem(<\/span>'lsqcurvefit'<\/span>, <\/span>...<\/span><\/span><\/div><\/div>
<\/span>'objective'<\/span>, model, <\/span>...<\/span><\/span><\/div><\/div>
<\/span>'xdata'<\/span>, t, <\/span>'ydata'<\/span>, c, <\/span>...<\/span><\/span><\/div><\/div>
<\/span>'x0'<\/span>,ones(1,6),<\/span>...<\/span><\/span><\/div><\/div>
<\/span>'lb'<\/span>, [-10 -10 -10 0 0 0 ],<\/span>...<\/span><\/span><\/div><\/div>
<\/span>'ub'<\/span>, [ 10 10 10 0.5 0.5 0.5], <\/span>...<\/span><\/span><\/div><\/div>
<\/span>'options'<\/span>,optimoptions(<\/span>'lsqcurvefit'<\/span>,<\/span>...<\/span><\/span><\/div><\/div>
<\/span>'OutputFcn'<\/span>, @curvefittingPlotIterates,<\/span>...<\/span><\/span><\/div><\/div>
<\/span>'Display'<\/span>,<\/span>'none'<\/span>))<\/span><\/span><\/div>
problem = struct with fields:<\/span><\/span><\/div>
objective: @(b,t)b(1)*exp(-b(4)*t)+b(2)*exp(-b(5)*t)+b(3)*exp(-b(6)*t)\r\n x0: [1 1 1 1 1 1]\r\n xdata: [3.9200 7.9300 11.8900 23.9000 47.8700 71.9100 93.8500 117.8400]\r\n ydata: [0.1630 0.6790 0.6790 0.3880 0.1830 0.1250 0.0860 0.0624]\r\n lb: [-10 -10 -10 0 0 0]\r\n ub: [10 10 10 0.5000 0.5000 0.5000]\r\n solver: 'lsqcurvefit'\r\n options: [1\u00d71 optim.options.Lsqcurvefit]<\/div><\/div><\/div><\/div><\/div><\/div>

Solve the Problem<\/span><\/h2>
First solve the problem directly once.<\/span><\/div>
b = lsqcurvefit(problem)<\/span><\/span><\/div>
<\/div><\/div>
b = <\/span>1\u00d76<\/span><\/div><\/div>
0.1842 0.1836 0.1841 0.0172 0.0171 0.0171\r\n<\/div>
<\/div>
<\/div><\/div><\/div><\/div><\/div><\/div><\/div>
You'll notice that the model does not do a stellar job fitting the data or even following the shape of the data. <\/span><\/div>

MultiStart<\/span><\/h2>
Let's see if we can do better by starting at a bunch of different points.<\/span><\/div>
ms = MultiStart;<\/span><\/span><\/div><\/div>
ms.Display = <\/span>'iter'<\/span>;<\/span><\/span><\/div><\/div>
rng <\/span>default<\/span><\/span><\/div><\/div>
figure<\/span><\/span><\/div><\/div>
tic<\/span><\/span><\/div><\/div>
[~,fval,exitflag,output,solutions] = run(ms, problem, 50)<\/span><\/span><\/div>
Run Local Local Local Local First-order\r\n Index exitflag f(x) # iter F-count optimality\r\n 1 3 0.222 8 63 0.0006396\r\n 2 3 0.000154 21 154 0.001864\r\n 3 3 0.009442 44 315 0.01989\r\n 4 3 1.462e-05 34 245 0.002586\r\n 5 3 1.454e-05 19 140 1.079e-05\r\n 6 3 1.475e-05 24 175 0.006883\r\n 7 3 0.009445 50 357 0.002266\r\n 8 3 1.495e-05 32 231 0.006853\r\n 9 3 1.466e-05 35 252 0.00478\r\n 10 3 0.00944 80 567 0.01042\r\n 11 3 1.471e-05 40 287 0.005472\r\n 12 3 1.566e-05 24 175 0.001576\r\n 13 3 0.009439 24 175 0.0005121\r\n 14 3 0.009451 41 294 0.02935\r\n 15 3 0.009493 26 189 0.004837\r\n 16 3 1.476e-05 40 287 0.006352\r\n 17 3 1.494e-05 40 287 0.008288\r\n 18 3 0.009446 62 441 0.01296\r\n 19 3 1.457e-05 22 161 0.001755\r\n 20 3 1.488e-05 53 378 0.007608\r\n 21 3 0.00944 37 266 0.006878\r\n 22 3 1.464e-05 24 175 0.003709\r\n 23 3 0.009449 43 308 0.02515\r\n 24 3 0.009447 47 336 0.007942\r\n 25 3 1.455e-05 23 168 0.001621\r\n 26 3 0.009442 32 231 0.01328\r\n 27 3 1.479e-05 40 287 0.004821\r\n 28 3 1.479e-05 18 133 0.006878\r\n 29 3 1.456e-05 72 511 0.0009721\r\n 30 3 1.455e-05 42 301 0.001122\r\n 31 3 0.009441 47 336 0.01537\r\n 32 3 0.009451 47 336 0.008942\r\n 33 3 0.0001729 14 105 0.0003276\r\n 34 3 0.009442 44 315 0.01062\r\n 35 3 0.0001751 21 154 7.71e-05\r\n 36 3 1.509e-05 26 189 0.009896\r\n 37 3 0.009458 39 280 0.02208\r\n 38 1 1.454e-05 24 175 7.815e-08\r\n 39 3 0.009441 60 427 0.0107\r\n 40 3 1.472e-05 34 245 0.002981\r\n 41 3 1.503e-05 22 161 0.00585\r\n 42 3 0.00952 15 112 0.008492\r\n 43 3 0.009439 21 154 0.000769\r\n 44 3 1.462e-05 64 455 0.0005576\r\n 45 3 0.009439 17 126 0.0001567\r\n 46 3 1.471e-05 30 217 0.001973\r\n 47 3 0.009444 38 273 0.02022\r\n 48 3 1.474e-05 24 175 0.004799\r\n 49 3 1.522e-05 42 301 0.008228<\/div><\/div>
<\/div><\/div>
50 3 0.009445 40 287 0.02166\r\n\r\nMultiStart completed the runs from all start points.\r\n\r\nAll 50 local solver runs converged with a positive local solver exit flag.<\/div><\/div>
fval = 1.4540e-05<\/div>
exitflag = 1<\/div>
output = struct with fields:<\/span><\/span><\/div>
funcCount: 12726\r\n localSolverTotal: 50\r\n localSolverSuccess: 50\r\n localSolverIncomplete: 0\r\n localSolverNoSolution: 0\r\n message: 'MultiStart completed the runs from all start points.\u21b5\u21b5All 50 local solver runs converged with a positive local solver exit flag.'<\/div><\/div><\/div>
solutions = <\/span>1\u00d750 object <\/span><\/div>
 <\/span><\/th>1<\/span><\/th>2<\/span><\/th>3<\/span><\/th>4<\/span><\/th>5<\/span><\/th>6<\/span><\/th>7<\/span><\/th>8<\/span><\/th>9<\/span><\/th>10<\/span><\/th>11<\/span><\/th>12<\/span><\/th>13<\/span><\/th>14<\/span><\/th>15<\/span><\/th>16<\/span><\/th>17<\/span><\/th>18<\/span><\/th>19<\/span><\/th>20<\/span><\/th>21<\/span><\/th>22<\/span><\/th>23<\/span><\/th>24<\/span><\/th>25<\/span><\/th>26<\/span><\/th>27<\/span><\/th>28<\/span><\/th>29<\/span><\/th>30<\/span><\/th><\/tr><\/thead>
1<\/span><\/th>1\u00d71 GlobalOptimSolution<\/span><\/td>1\u00d71 GlobalOptimSolution<\/span><\/td>1\u00d71 GlobalOptimSolution<\/span><\/td>1\u00d71 GlobalOptimSolution<\/span><\/td>1\u00d71 GlobalOptimSolution<\/span><\/td>1\u00d71 GlobalOptimSolution<\/span><\/td>1\u00d71 GlobalOptimSolution<\/span><\/td>1\u00d71 GlobalOptimSolution<\/span><\/td>1\u00d71 GlobalOptimSolution<\/span><\/td>1\u00d71 GlobalOptimSolution<\/span><\/td>1\u00d71 GlobalOptimSolution<\/span><\/td>1\u00d71 GlobalOptimSolution<\/span><\/td>1\u00d71 GlobalOptimSolution<\/span><\/td>1\u00d71 GlobalOptimSolution<\/span><\/td>1\u00d71 GlobalOptimSolution<\/span><\/td>1\u00d71 GlobalOptimSolution<\/span><\/td>1\u00d71 GlobalOptimSolution<\/span><\/td>1\u00d71 GlobalOptimSolution<\/span><\/td>1\u00d71 GlobalOptimSolution<\/span><\/td>1\u00d71 GlobalOptimSolution<\/span><\/td>1\u00d71 GlobalOptimSolution<\/span><\/td>1\u00d71 GlobalOptimSolution<\/span><\/td>1\u00d71 GlobalOptimSolution<\/span><\/td>1\u00d71 GlobalOptimSolution<\/span><\/td>1\u00d71 GlobalOptimSolution<\/span><\/td>1\u00d71 GlobalOptimSolution<\/span><\/td>1\u00d71 GlobalOptimSolution<\/span><\/td>1\u00d71 GlobalOptimSolution<\/span><\/td>1\u00d71 GlobalOptimSolution<\/span><\/td>1\u00d71 GlobalOptimSolution<\/span><\/td><\/tr><\/tbody><\/table><\/div><\/div><\/div><\/div><\/div>
serialTime = toc;<\/span><\/span><\/div><\/div><\/div>

Visualize the Best Solution<\/span><\/h2>
The 50th solution, which is what is plotted above, is not necessarily the best one. Luckily for us, <\/span>MultiStart<\/span><\/a> orders the solutions from best to worst. So we need only look at the first one.<\/span><\/div>
curvefittingPlotIterates(solutions)<\/span><\/span><\/div>
<\/div><\/div><\/div><\/div><\/div>
You can see now that the 50th was not the best solution as the mean squared error on this final one displayed is over a factor of 10 better.<\/span><\/div>

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

Calculate Speedup<\/span><\/h2>
Speed up may not be evident until second run due to pool start up time. Since I started mine earlier, I get to see decent speed up.<\/span><\/div>
speedup = serialTime\/parallelTime<\/span><\/span><\/div>
speedup = 2.5014<\/div><\/div><\/div><\/div>

Do You Have Problems Where the Solution is Sensitive to the Starting Point?<\/span><\/h2>
Tell us<\/span><\/a> about your exploration of the solution space to your problem. If the solution is sensitive to where you start, you might consider using <\/span>MultiStart<\/span> and other techniques from the Global Optimization Toolbox.<\/span><\/div>
Copyright 2021 The MathWorks, Inc.<\/span><\/div>

Appendix<\/span><\/h2>
Here's the code for plotting the iterates.<\/span><\/div>
dbtype <\/span>curvefittingPlotIterates<\/span><\/span><\/div>
1 function stop = curvefittingPlotIterates(x,optimValues,state)\r\n2 % Output function that plots the iterates of the optimization algorithm.\r\n3 \r\n4 % Copyright 2010 The MathWorks, Inc.\r\n5 \r\n6 persistent x0 r;\r\n7 if nargin == 1\r\n8 showPlot(x(1).X,x(1).X0{:},x(1).Fval)\r\n9 else\r\n10 switch state\r\n11 case 'init' % store initial point for later use\r\n12 x0 = x;\r\n13 case 'done'\r\n14 if ~(optimValues.iteration == 0)\r\n15 % After optimization, display solution in plot title\r\n16 r = optimValues.resnorm;\r\n17 showPlot(x,x0,r)\r\n18 end\r\n19 end\r\n20 end\r\n21 if nargout > 0\r\n22 stop = false;\r\n23 clear function\r\n24 end\r\n25 end\r\n26 \r\n27 function showPlot(b,b0,r)\r\n28 f = @(b,x) b(1)*exp(-b(4).*x) + b(2).*exp(-b(5).*x) +...\r\n29 b(3).*exp(-b(6).*x);\r\n30 \r\n31 persistent h ha\r\n32 if isempty(h) || ~isvalid(h)\r\n33 x = [ 3.92, 7.93, 11.89, 23.90, 47.87, 71.91, 93.85, 117.84 ];\r\n34 y = [ 0.163, 0.679, 0.679, 0.388, 0.183, 0.125, 0.086, 0.0624 ];\r\n35 plot(x,y,'o');\r\n36 xlabel('t')\r\n37 ylabel('c')\r\n38 title('c=b_1e^{-b_4t}+b_2e^{-b_5t}+b_3e^{-b_6t}')\r\n39 axis([0 120 0 0.8]); \r\n40 h = line(3:120,f(b,3:120),'Color','r','Tag','PlotIterates');\r\n41 \r\n42 else\r\n43 set(h,'YData',f(b,get(h,'XData')));\r\n44 end\r\n45 s = sprintf('Starting Value Fitted Value\\n\\n');\r\n46 \r\n47 for i = 1:length(b)\r\n48 s = [s, sprintf('b(%d): % 2.4f b(%d): % 2.4f\\n',i,b0(i),i,b(i))];\r\n49 end\r\n50 s = [s,sprintf('\\nMSE = %2.4e',r)];\r\n51 \r\n52 if isempty(ha) || ~isvalid(ha)\r\n53 % Create textbox\r\n54 ha = annotation(gcf,'textbox',...\r\n55 [0.5 0.5 0.31 0.32],...\r\n56 'String',s,...\r\n57 'FitBoxToText','on',...\r\n58 'Tag','CoeffDisplay');\r\n59 end\r\n60 ha.String = s;\r\n61 drawnow\r\n62 \r\n63 end<\/div><\/div><\/div><\/div><\/div>\r\n<\/div>