A few weeks ago when Will picked Wing Designer 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 hand, I got a score of 30706.9.
Fortunately with the way John wrote this, I could tie into the scoring engine by emulating the output of user selections from the user interface.
Looking at the tunable parameters there are quite a few and many of them are integer constrained. Since the actual engine for this app is a black box to me, I am assuming that the score function is highly nonlinear with lots of local minima. Thus I need to take a global approach to this. The first function that jumps out is ga since this allows integer constraints. However, genetic algorithms don't always hold up well in as many dimensions as we have.
So instead I turned to patternsearch. Patternsearch is a direct search method that is often recommended over ga by Alan, who's a genius. But since patternsearch does not support integer constraints, I need to evaluate it on every combination of integers and then pick the best from these results.
Wing Designer's user interface returns a geometry structure that I needed to emulate. Here's an example of part of the structure:
I wrote a function that takes in some of the various parameters that would be selectable in the user interface to create this structure.
Next, I needed to make some design decisions regarding my wing since trying every integer combination is just not feasible in the amount of time allotted. I spoke with one of our aerospace experts and he provided a few suggestions:
- Wings should be long and narrow.
- Don't break the sound barrier.
- Probably use a jet.
- Fly between 35k and 50k.
- Some of the parameters like friction would be known before hand so don't include them in the optimization.
I made an additional few:
- Root and Tip should have the same airfoil.
- I will use 10 spanwise panels and 6 chordwise panels. More than this dramatically slows the computation.
Since patternsearch has to run many times, once for each of the airfoils, and the function is called many times inside of patternsearch, I have to parallelize it. I could set the 'UseParallel' flag in psoptimset, but this would only parallelize individual patternsearch polls. Ideally, the calls to patternsearch should be parallelized, so instead we'll use a parallel for-loop (parfor) over the airfoils.
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 a cluster running MATLAB Distributed Computing Server on Amazon EC2. I fired up a cluster with 36 workers which makes it much more feasible. I chose 36 workers: one for each airfoil, one to run the job itself, and an extra because it likes even numbers.
This is the body of OptimizeWing.m, the function that calls patternsearch.
% Lower and upper bounds lb = [50 10 5 0 0 0 0 300 35000 0 0 5000 0]; ub = [200 40 10 5 5 5 5 660 50000 7 15 80000 100]; % x0, half way between bounds x0 = lb + ((ub-lb)./2); % Options options = psoptimset('TolMesh',0.1,'MaxIter',100); nfoils = 34; % 34 airfoils fval = zeros(nfoils,1); % fval at xopt xopt = zeros(nfoils,numel(x0)); % optimal x % Run the parfor loop parfor ii = 1:nfoils % Wrap the call in a try block in case some air foils cause failure try [xopt(ii,:), fval(ii)] = patternsearch(@(x)wingdesignerScore(x,ii),x0,,,,,lb,ub,,options) catch % iith airfoil failed xopt(ii,:) = nan; fval(ii) = nan; end end
First, I ran patternsearch once on one airfoil with very loose tolerances just to make sure it worked as expected. It did! Then it was time to submit the main function to the cluster, using batch.
And wait... Well actually, shut off my computer and go home for the evening. Since the job is running on EC2 somewhere in Virginia, there is no dependency on my machine.
After retrieving the completed job in the morning, the best score is found by taking the minimum of the scores vector, fval and putting the corresponding xopt parameters into the App.
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?
Do you have a challenging optimization problem or some really computationally expensive work where cluster or cloud computing could help?
Let us know about it in the comments!
Get the MATLAB code
Published with MATLAB® R2014a
2 CommentsOldest to Newest
One issue is that the aircraft actually is supersonic – the airspeed bound is for speed of sound at sea level, which is higher than that at 35,000 feet. A wing like that probably wouldn’t fare too well crossing the sound barrier. Also, the aircraft’s empty weight doesn’t seem to couple into anything else, so naturally, the optimizer settled onto the lowest possible value, producing a fuel fraction of 0.6%. For reference, the Rutan Voyager, which flew around the world in 1986, had a fuel fraction of about 10%, and that thing was about as fragile as could be.
We just finished up a student design project for a long-endurance UAV, so I had a good chuckle when I saw this pick show up a few days after submission. We ended up with a decent little sizing algorithm, but since no one knew much of anything about optimization, we resorted to the brute force method of computing the entire design domain. (We were looking for minimum airframe cost for a given cruise speed, endurance, and wingspan, so it wasn’t prohibitive.) With another couple of weeks, we probably could have plugged the routine into one of the optimization functions, but there just wasn’t quite enough time in the semester for that extra bit of fun.
I saw that the value in the app did show greater than one mach 1. I guess I could have added a nonlinear constraint to insure this but that would’ve required many more function evaluations. I didn’t look into the Wing Designer code enough to see if properly punished breaking the sound barrier.
It did have a few additional checks inside of the app for parameters such as payload (what allows you to have a negative score) so I am guessing that adding a few more could keep the optimizer in check without even needing nonlinear constraints.