{"id":5341,"date":"2014-05-30T09:00:32","date_gmt":"2014-05-30T13:00:32","guid":{"rendered":"https:\/\/blogs.mathworks.com\/pick\/?p=5341"},"modified":"2018-03-15T16:31:45","modified_gmt":"2018-03-15T20:31:45","slug":"wing-designer-part-ii-optimize","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/pick\/2014\/05\/30\/wing-designer-part-ii-optimize\/","title":{"rendered":"Wing Designer Part II- Optimize!"},"content":{"rendered":"<div xmlns:mwsh=\"https:\/\/www.mathworks.com\/namespace\/mcode\/v1\/syntaxhighlight.dtd\" class=\"content\">\r\n   <introduction>\r\n      <p><a href=\"https:\/\/www.mathworks.com\/matlabcentral\/answers\/contributors\/3208495\">Sean<\/a>'s pick this week is also <a href=\"https:\/\/www.mathworks.com\/matlabcentral\/fileexchange\/15442-wing-designer\">Wing Designer<\/a> by <a href=\"https:\/\/www.mathworks.com\/matlabcentral\/fileexchange\/authors\/28606\">John Rogers<\/a>.\r\n      <\/p>\r\n   <\/introduction>\r\n   <h3>Contents<\/h3>\r\n   <div>\r\n      <ul>\r\n         <li><a href=\"#1\">Background<\/a><\/li>\r\n         <li><a href=\"#2\">Picking an Optimizer<\/a><\/li>\r\n         <li><a href=\"#3\">Implementation<\/a><\/li>\r\n         <li><a href=\"#5\">Making this Feasible<\/a><\/li>\r\n         <li><a href=\"#8\">Results<\/a><\/li>\r\n         <li><a href=\"#9\">Comments<\/a><\/li>\r\n      <\/ul>\r\n   <\/div>\r\n   <h3>Background<a name=\"1\"><\/a><\/h3>\r\n   <p>A few weeks ago when <a href=\"https:\/\/blogs.mathworks.com\/pick\/2014\/05\/16\/wing-designer\/\">Will picked Wing Designer<\/a> I immediately had to go download it and try to beat his high score of 29252.  After a bit of twiddling I was able to!  By\r\n      hand, I got a score of 30706.9.\r\n   <\/p>\r\n   <p>But then I wondered what kind of score I could get with an optimizer from the <a href=\"https:\/\/www.mathworks.com\/products\/optimization\/\">Optimization Toolbox<\/a> or the <a href=\"https:\/\/www.mathworks.com\/products\/global-optimization\/\">Global Optimization Toolbox<\/a>?\r\n   <\/p>\r\n   <p>Fortunately with the way John wrote this, I could tie into the scoring engine by emulating the output of user selections from\r\n      the user interface.\r\n   <\/p>\r\n   <h3>Picking an Optimizer<a name=\"2\"><\/a><\/h3>\r\n   <p>Looking at the tunable parameters there are quite a few and many of them are integer constrained.  Since the actual engine\r\n      for this app is a black box to me, I am <i>assuming<\/i> that the score function is highly nonlinear with lots of local minima.  Thus I need to take a global approach to this.  The\r\n      first function that jumps out is <a href=\"\"><tt>ga<\/tt><\/a> since this allows integer constraints.  However, genetic algorithms don't always hold up well in as many dimensions as we\r\n      have.\r\n   <\/p>\r\n   <p>So instead I turned to <a href=\"\"><tt>patternsearch<\/tt><\/a>.  Patternsearch is a direct search method that is often <a href=\"https:\/\/www.mathworks.com\/matlabcentral\/answers\/75753-suggestion-for-using-global-optimization#answer_85804\">recommended<\/a> over <tt>ga<\/tt> by <a href=\"https:\/\/www.mathworks.com\/matlabcentral\/answers\/contributors\/1033975-alan-weiss\">Alan<\/a>, who's a genius.  But since <tt>patternsearch<\/tt> does not support integer constraints, I need to evaluate it on every combination of integers and then pick the best from\r\n      these results.\r\n   <\/p>\r\n   <h3>Implementation<a name=\"3\"><\/a><\/h3>\r\n   <p>Wing Designer's user interface returns a geometry structure that I needed to emulate.  Here's an example of part of the structure:<\/p>\r\n   <p><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/pick\/Sean\/mainwd\/geostruct.PNG\"> <\/p>\r\n   <p>I wrote a function that takes in some of the various parameters that would be selectable in the user interface to create this\r\n      structure.\r\n   <\/p>\r\n   <p>Next, I needed to make some design decisions regarding my wing since trying every integer combination is just not feasible\r\n      in the amount of time allotted.  I spoke with one of our aerospace experts and he provided a few suggestions:\r\n   <\/p>\r\n   <div>\r\n      <ul>\r\n         <li>Wings should be long and narrow.<\/li>\r\n         <li>Don't break the sound barrier.<\/li>\r\n         <li>Probably use a jet.<\/li>\r\n         <li>Fly between 35k and 50k.<\/li>\r\n         <li>Some of the parameters like friction would be known before hand so   don't include them in the optimization.<\/li>\r\n      <\/ul>\r\n   <\/div>\r\n   <p>I made an additional few:<\/p>\r\n   <div>\r\n      <ul>\r\n         <li>Root and Tip should have the same airfoil.<\/li>\r\n         <li>I will use 10 spanwise panels and 6 chordwise panels.  More than this dramatically slows the computation.<\/li>\r\n      <\/ul>\r\n   <\/div>\r\n   <h3>Making this Feasible<a name=\"5\"><\/a><\/h3>\r\n   <p>Since <tt>patternsearch<\/tt> has to run many times, once for each of the airfoils, and the function is called many times inside of <tt>patternsearch<\/tt>, I have to parallelize it.  I could set the <i>'UseParallel'<\/i> flag in <tt>psoptimset<\/tt>, but this would only parallelize individual <tt>patternsearch<\/tt> polls. Ideally, the calls to <tt>patternsearch<\/tt> should be parallelized, so instead we'll use a <a href=\"\">parallel for-loop (parfor)<\/a> over the airfoils.\r\n   <\/p>\r\n   <p>My laptop only has two cores, not enough to make this happen in the 12 hours before the blog post is due.  So instead, I used\r\n      a cluster running <a href=\"https:\/\/www.mathworks.com\/products\/distriben\/\">MATLAB Distributed Computing Server<\/a> on <a href=\"https:\/\/www.mathworks.com\/products\/parallel-computing\/parallel-computing-on-the-cloud\/distriben-ec2.html\">Amazon EC2<\/a>.  I fired up a cluster with 36 workers which makes it much more feasible.  I chose 36 workers: one for each airfoil, one\r\n      to run the job itself, and an extra because it likes even numbers.\r\n   <\/p>\r\n   <p><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/pick\/Sean\/mainwd\/ec2cluster.PNG\"> <\/p>\r\n   <p>This is the body of <tt>OptimizeWing.m<\/tt>, the function that calls <tt>patternsearch<\/tt>.\r\n   <\/p>\r\n   <p><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/pick\/Sean\/mainwd\/xlabels.PNG\"> <\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\"><span style=\"color: #228B22\">% Lower and upper bounds<\/span>\r\nlb = [50 10 5 0 0 0 0 300 35000 0 0 5000 0];\r\nub = [200 40 10 5 5 5 5 660 50000 7 15 80000 100];\r\n\r\n<span style=\"color: #228B22\">% x0, half way between bounds<\/span>\r\nx0 = lb + ((ub-lb).\/2);\r\n\r\n<span style=\"color: #228B22\">% Options<\/span>\r\noptions = psoptimset(<span style=\"color: #A020F0\">'TolMesh'<\/span>,0.1,<span style=\"color: #A020F0\">'MaxIter'<\/span>,100);\r\n\r\nnfoils = 34; <span style=\"color: #228B22\">% 34 airfoils<\/span>\r\nfval = zeros(nfoils,1); <span style=\"color: #228B22\">% fval at xopt<\/span>\r\nxopt = zeros(nfoils,numel(x0)); <span style=\"color: #228B22\">% optimal x<\/span>\r\n\r\n<span style=\"color: #228B22\">% Run the parfor loop<\/span>\r\n<span style=\"color: #0000FF\">parfor<\/span> ii = 1:nfoils\r\n    <span style=\"color: #228B22\">% Wrap the call in a try block in case some air foils cause failure<\/span>\r\n    <span style=\"color: #0000FF\">try<\/span>\r\n        [xopt(ii,:), fval(ii)] = patternsearch(@(x)wingdesignerScore(x,ii),x0,[],[],[],[],lb,ub,[],options)\r\n    <span style=\"color: #0000FF\">catch<\/span>\r\n        <span style=\"color: #228B22\">% iith airfoil failed<\/span>\r\n        xopt(ii,:) = nan;\r\n        fval(ii) = nan;\r\n    <span style=\"color: #0000FF\">end<\/span>\r\n<span style=\"color: #0000FF\">end<\/span><\/pre><p>First, I ran <tt>patternsearch<\/tt> once on one airfoil with very loose tolerances just to make sure it worked as expected.  It did!  Then it was time to submit\r\n      the main function to the cluster, using <a href=\"\"><tt>batch<\/tt><\/a>.\r\n   <\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">batch(<span style=\"color: #A020F0\">'OptimizeWing.m'<\/span>,2,<span style=\"color: #A020F0\">'Pool'<\/span>,34,<span style=\"color: #A020F0\">'AttachedFiles'<\/span>,files)<\/pre><p>And wait... Well actually, shut off my computer and go home for the evening. Since the job is running on EC2 somewhere in\r\n      Virginia, there is no dependency on my machine.\r\n   <\/p>\r\n   <h3>Results<a name=\"8\"><\/a><\/h3>\r\n   <p>After retrieving the completed job in the morning, the best score is found by taking the minimum of the scores vector, <i>fval<\/i> and putting the corresponding <i>xopt<\/i> parameters into the App.\r\n   <\/p>\r\n   <p><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/pick\/Sean\/mainwd\/AppFinal.PNG\"><br><br>That's about 153000x better than we did by hand!  The optimized wing is large, symmetric, and mostly filled with fuel.  I'd be curious to hear thoughts from someone in the industry on the feasibility of a wing like this?<\/p>\r\n   <h3>Comments<a name=\"9\"><\/a><\/h3>\r\n   <p>Do you have a challenging optimization problem or some really computationally expensive work where cluster or cloud computing\r\n      could help?\r\n   <\/p>\r\n   <p>Let us know about it in the <a href=\"https:\/\/blogs.mathworks.com\/pick\/?p=5341#respond\">comments<\/a>!\r\n   <\/p><script language=\"JavaScript\">\r\n<!--\r\n\r\n    function grabCode_54ced8c0cb4e48dfb6f2e048b8e955a8() {\r\n        \/\/ Remember the title so we can use it in the new page\r\n        title = document.title;\r\n\r\n        \/\/ Break up these strings so that their presence\r\n        \/\/ in the Javascript doesn't mess up the search for\r\n        \/\/ the MATLAB code.\r\n        t1='54ced8c0cb4e48dfb6f2e048b8e955a8 ' + '##### ' + 'SOURCE BEGIN' + ' #####';\r\n        t2='##### ' + 'SOURCE END' + ' #####' + ' 54ced8c0cb4e48dfb6f2e048b8e955a8';\r\n    \r\n        b=document.getElementsByTagName('body')[0];\r\n        i1=b.innerHTML.indexOf(t1)+t1.length;\r\n        i2=b.innerHTML.indexOf(t2);\r\n \r\n        code_string = b.innerHTML.substring(i1, i2);\r\n        code_string = code_string.replace(\/REPLACE_WITH_DASH_DASH\/g,'--');\r\n\r\n        \/\/ Use \/x3C\/g instead of the less-than character to avoid errors \r\n        \/\/ in the XML parser.\r\n        \/\/ Use '\\x26#60;' instead of '<' so that the XML parser\r\n        \/\/ doesn't go ahead and substitute the less-than character. \r\n        code_string = code_string.replace(\/\\x3C\/g, '\\x26#60;');\r\n\r\n        author = 'Sean de Wolski';\r\n        copyright = 'Copyright 2014 The MathWorks, Inc.';\r\n\r\n        w = window.open();\r\n        d = w.document;\r\n        d.write('<pre>\\n');\r\n        d.write(code_string);\r\n\r\n        \/\/ Add author and copyright lines at the bottom if specified.\r\n        if ((author.length > 0) || (copyright.length > 0)) {\r\n            d.writeln('');\r\n            d.writeln('%%');\r\n            if (author.length > 0) {\r\n                d.writeln('% _' + author + '_');\r\n            }\r\n            if (copyright.length > 0) {\r\n                d.writeln('% _' + copyright + '_');\r\n            }\r\n        }\r\n\r\n        d.write('<\/pre>\\n');\r\n      \r\n      d.title = title + ' (MATLAB code)';\r\n      d.close();\r\n      }   \r\n      \r\n-->\r\n<\/script><p style=\"text-align: right; font-size: xx-small; font-weight:lighter;   font-style: italic; color: gray\"><br><a href=\"javascript:grabCode_54ced8c0cb4e48dfb6f2e048b8e955a8()\"><span style=\"font-size: x-small;        font-style: italic;\">Get \r\n            the MATLAB code \r\n            <noscript>(requires JavaScript)<\/noscript><\/span><\/a><br><br>\r\n      Published with MATLAB&reg; R2014a<br><\/p>\r\n<\/div>\r\n<!--\r\n54ced8c0cb4e48dfb6f2e048b8e955a8 ##### SOURCE BEGIN #####\r\n%% Wing Designer - Part 2\r\n%\r\n% <https:\/\/www.mathworks.com\/matlabcentral\/answers\/contributors\/3208495\r\n% Sean>'s pick this week is also\r\n% <https:\/\/www.mathworks.com\/matlabcentral\/fileexchange\/15442-wing-designer Wing Designer> by\r\n% <https:\/\/www.mathworks.com\/matlabcentral\/fileexchange\/authors\/28606 John Rogers>.\r\n% \r\n% \r\n\r\n%% Background\r\n% A few weeks ago when\r\n% <https:\/\/blogs.mathworks.com\/pick\/2014\/05\/16\/wing-designer\/ Will picked\r\n% Wing Designer> I immediately had to go download it and try to beat his\r\n% high score of 29252.  After a bit of twiddling I was able to!  By hand, I\r\n% got a score of 30706.9.\r\n%\r\n% But then I wondered what kind of score I could get with an optimizer from\r\n% the <https:\/\/www.mathworks.com\/products\/optimization\/ Optimization\r\n% Toolbox> or the <https:\/\/www.mathworks.com\/products\/global-optimization\/\r\n% Global Optimization Toolbox>?\r\n%\r\n% Fortunately with the way John wrote this, I could tie into the scoring\r\n% engine by emulating the output of user selections from the user interface.\r\n%\r\n\r\n%% Picking an Optimizer\r\n% Looking at the tunable parameters there are quite a few and many of them\r\n% are integer constrained.  Since the actual engine for this app is a black\r\n% box to me, I am _assuming_ that the score function is highly nonlinear\r\n% with lots of local minima.  Thus I need to take a global approach to\r\n% this.  The first function that jumps out is\r\n% < |ga|> since\r\n% this allows integer constraints.  However, genetic algorithms don't\r\n% always hold up well in as many dimensions as we have.\r\n%\r\n% So instead I turned to\r\n% <\r\n% |patternsearch|>.  Patternsearch is a direct search method that is often\r\n% <https:\/\/www.mathworks.com\/matlabcentral\/answers\/75753-suggestion-for-using-global-optimization#answer_85804\r\n% recommended> over |ga| by\r\n% <https:\/\/www.mathworks.com\/matlabcentral\/answers\/contributors\/1033975-alan-weiss\r\n% Alan>, who's a genius.  But since |patternsearch| does not support\r\n% integer constraints, I need to evaluate it on every combination of\r\n% integers and then pick the best from these results.\r\n%\r\n\r\n%% Implementation\r\n% Wing Designer's user interface returns a geometry structure that I needed\r\n% to emulate.  Here's an example of part of the structure:\r\n%\r\n% <<geostruct.PNG>>\r\n\r\n\r\n%% \r\n% I wrote a function that takes in some of the various parameters that\r\n% would be selectable in the user interface to create this structure.\r\n%\r\n% Next, I needed to make some design decisions regarding my wing since\r\n% trying every integer combination is just not feasible in the amount of\r\n% time allotted.  I spoke with one of our aerospace experts and he provided\r\n% a few suggestions:\r\n%\r\n% * Wings should be long and narrow.\r\n% * Don't break the sound barrier.\r\n% * Probably use a jet.\r\n% * Fly between 35k and 50k.\r\n% * Some of the parameters like friction would be known before hand so\r\n%   don't include them in the optimization.\r\n% \r\n% I made an additional few:\r\n% \r\n% * Root and Tip should have the same airfoil.\r\n% * I will use 10 spanwise panels and 6 chordwise panels.  More than this\r\n% dramatically slows the computation.\r\n%\r\n\r\n%% Making this Feasible\r\n% \r\n% Since |patternsearch| has to run many times, once for each of the\r\n% airfoils, and the function is called many times inside of\r\n% |patternsearch|, I have to parallelize it.  I could set the\r\n% _'UseParallel'_ flag in |psoptimset|, but this would only parallelize\r\n% individual |patternsearch| polls. Ideally, the calls to |patternsearch|\r\n% should be parallelized, so instead we'll use a\r\n% <\r\n% parallel for-loop (parfor)> over the airfoils.\r\n%\r\n% My laptop only has two cores, not enough to make this happen in the 12\r\n% hours before the blog post is due.  So instead, I used a cluster running\r\n% <https:\/\/www.mathworks.com\/products\/distriben\/ MATLAB Distributed\r\n% Computing Server> on <https:\/\/www.mathworks.com\/products\/parallel-computing\/parallel-computing-on-the-cloud\/distriben-ec2.html\r\n% Amazon EC2>.  I fired up a cluster with 36 workers which makes it much\r\n% more feasible.  I chose 36 workers: one for each airfoil, one to run the\r\n% job itself, and an extra because it likes even numbers.\r\n%\r\n% <<ec2cluster.PNG>>\r\n%\r\n% This is the body of |OptimizeWing.m|, the function that calls\r\n% |patternsearch|.\r\n%\r\n% <<xlabels.PNG>>\r\n\r\n% Lower and upper bounds\r\nlb = [50 10 5 0 0 0 0 300 35000 0 0 5000 0];\r\nub = [200 40 10 5 5 5 5 660 50000 7 15 80000 100];    \r\n  \r\n% x0, half way between bounds\r\nx0 = lb + ((ub-lb).\/2);\r\n\r\n% Options\r\noptions = psoptimset('TolMesh',0.1,'MaxIter',100);\r\n\r\nnfoils = 34; % 34 airfoils\r\nfval = zeros(nfoils,1); % fval at xopt\r\nxopt = zeros(nfoils,numel(x0)); % optimal x\r\n\r\n% Run the parfor loop\r\nparfor ii = 1:nfoils\r\n    % Wrap the call in a try block in case some air foils cause failure\r\n    try\r\n        [xopt(ii,:), fval(ii)] = patternsearch(@(x)wingdesignerScore(x,ii),x0,[],[],[],[],lb,ub,[],options)\r\n    catch\r\n        % iith airfoil failed\r\n        xopt(ii,:) = nan;\r\n        fval(ii) = nan;\r\n    end\r\nend\r\n\r\n\r\n%% \r\n% First, I ran |patternsearch| once on one airfoil with very loose\r\n% tolerances just to make sure it worked as expected.  It did!  Then it was\r\n% time to submit the main function to the cluster, using\r\n% <\r\n% |batch|>.\r\n%\r\n\r\nbatch('OptimizeWing.m',2,'Pool',34,'AttachedFiles',files)\r\n\r\n\r\n%% \r\n% And wait... Well actually, shut off my computer and go home for the\r\n% evening. Since the job is running on EC2 somewhere in Virginia, there is\r\n% no dependency on my machine.\r\n%\r\n\r\n\r\n%% Results\r\n% After retrieving the completed job in the morning, the best score is\r\n% found by taking the minimum of the scores vector, _fval_ and putting the\r\n% corresponding _xopt_ parameters into the App.\r\n%\r\n% <<AppFinal.PNG>>\r\n%\r\n\r\n%% Comments\r\n% \r\n% Do you have a challenging optimization problem or some really\r\n% computationally expensive work where cluster or cloud computing could\r\n% help?\r\n%\r\n% Let us know about it in the <https:\/\/blogs.mathworks.com\/pick\/?p=5341#respond\r\n% comments>! \r\n%\r\n\r\n\r\n##### SOURCE END ##### 54ced8c0cb4e48dfb6f2e048b8e955a8\r\n-->","protected":false},"excerpt":{"rendered":"<div class=\"overview-image\"><img decoding=\"async\"  class=\"img-responsive\" src=\"https:\/\/blogs.mathworks.com\/images\/pick\/Sean\/mainwd\/geostruct.PNG\" onError=\"this.style.display ='none';\" \/><\/div><p>\r\n   \r\n      Sean's pick this week is also Wing Designer by John Rogers.\r\n      \r\n   \r\n   Contents\r\n   \r\n      \r\n        ... <a class=\"read-more\" href=\"https:\/\/blogs.mathworks.com\/pick\/2014\/05\/30\/wing-designer-part-ii-optimize\/\">read more >><\/a><\/p>","protected":false},"author":87,"featured_media":0,"comment_status":"open","ping_status":"closed","sticky":false,"template":"","format":"standard","meta":[],"categories":[16],"tags":[],"_links":{"self":[{"href":"https:\/\/blogs.mathworks.com\/pick\/wp-json\/wp\/v2\/posts\/5341"}],"collection":[{"href":"https:\/\/blogs.mathworks.com\/pick\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/blogs.mathworks.com\/pick\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/blogs.mathworks.com\/pick\/wp-json\/wp\/v2\/users\/87"}],"replies":[{"embeddable":true,"href":"https:\/\/blogs.mathworks.com\/pick\/wp-json\/wp\/v2\/comments?post=5341"}],"version-history":[{"count":11,"href":"https:\/\/blogs.mathworks.com\/pick\/wp-json\/wp\/v2\/posts\/5341\/revisions"}],"predecessor-version":[{"id":9564,"href":"https:\/\/blogs.mathworks.com\/pick\/wp-json\/wp\/v2\/posts\/5341\/revisions\/9564"}],"wp:attachment":[{"href":"https:\/\/blogs.mathworks.com\/pick\/wp-json\/wp\/v2\/media?parent=5341"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/pick\/wp-json\/wp\/v2\/categories?post=5341"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/pick\/wp-json\/wp\/v2\/tags?post=5341"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}