diff --git a/README.md b/README.md new file mode 100644 index 0000000..920dcac --- /dev/null +++ b/README.md @@ -0,0 +1,15 @@ +# PlanetaryGearbox + +To run the planetary gearbox optimization file, the required files to download are: + +constraint.m +evalFitness.m +explainX.m +fval.m + +GA550.m +goptions.m +runAllOptimizations.m +runOneOptimization.m + +When all the files are in the same folder in MATLAB, open runAllOptimizations, set the desired minimum constraints for the sun and ring gears, and run the script. diff --git a/constraint.m b/constraint.m index 0a94525..4c2e1cf 100644 --- a/constraint.m +++ b/constraint.m @@ -1,57 +1,84 @@ -function val = constraint(x, bearing_d, sun_extrusion_d) -%% Planetary Gearbox GA Constraint Function -% Elliot Stockwell 1/12/21 +function val = constraint(x, bearing_diameter, planet_shaft_diameter, sun_shaft_diameter, maxX) +%% Gear Ratio Constraint Function % -% x : array of input variables (planet_d1, planet_d2, sun_d1, dp) -% bearing_d : inner diameter of the bearing +% X is formatted +% [ +% front_ring_teeth / 3 +% front_planet_s2_teeth +% front_sun_teeth / 3 +% rear_planet_s2_teeth +% rear_sun_teeth / 3 +% diametral_pitch +% ] -pm = 1E2; % penalty multiplier +%% GA Constants +PENALTY = 1e2; +val = 0; -planet_d1 = x(1); -planet_d2 = x(2); -sun_d1 = x(3); -dp = 12 + 4 * x(4); -ring_d2 = planet_d1 + planet_d2 + sun_d1; +%% Interpret x +front_ring_teeth = 3 * x(1); +front_planet_s2_teeth = x(2); +front_sun_teeth = 3 * x(3); +rear_planet_s2_teeth = x(4); +rear_sun_teeth = 3 * x(5); +dp = 12 + 4 * x(6); -g = zeros(1, 8); -mod_tolerance = .01; +%% Calculate derived variables +rear_ring_teeth = front_ring_teeth; -%% Make sure teeth are integers -planet_1_mod = mod(dp * planet_d1, 1); -planet_2_mod = mod(dp * planet_d2, 1); -sun_mod = mod(dp * sun_d1, 3); -ring_mod = mod(dp * ring_d2, 3); +front_planet_s1_teeth = front_ring_teeth - ... + front_planet_s2_teeth - front_sun_teeth; -if sun_mod < 3 - mod_tolerance && sun_mod > mod_tolerance - g(1) = pm; -end +rear_planet_s1_teeth = rear_ring_teeth - ... + rear_planet_s2_teeth - rear_sun_teeth; -if ring_mod < 3 - mod_tolerance && ring_mod > mod_tolerance - g(2) = pm; -end -if planet_1_mod < 1 - mod_tolerance && planet_1_mod > mod_tolerance - g(3) = pm; +%% Calculate diameters +front_ring_diameter = front_ring_teeth / dp; +front_planet_s1_diameter = front_planet_s1_teeth / dp; +front_planet_s2_diameter = front_planet_s2_teeth / dp; +front_sun_diameter = front_sun_teeth / dp; + +rear_ring_diameter = rear_ring_teeth / dp; +rear_planet_s1_diameter = rear_planet_s1_teeth / dp; +rear_planet_s2_diameter = rear_planet_s2_teeth / dp; +rear_sun_diameter = rear_sun_teeth / dp; + + +%% Calculate side constraints +% Variables +comp_planet_diameter = (front_sun_diameter + front_planet_s1_diameter) / (rear_sun_diameter + rear_planet_s1_diameter) - 1; + +% Don't exceed bearing bounds +if 2 * front_planet_s1_diameter + front_sun_diameter > bearing_diameter + val = val + PENALTY; end -if planet_2_mod < 1 - mod_tolerance && planet_2_mod > mod_tolerance - g(4) = pm; +if 2 * rear_planet_s1_diameter + rear_sun_diameter > bearing_diameter + val = val + PENALTY; end -%% Side constraints -if dp > 24 - g(5) = pm; +if comp_planet_diameter < -0.01 || comp_planet_diameter > 0.01 + val = val + PENALTY; end -g(6) = max(0, (sun_d1 + 2 * planet_d1) / bearing_d - 1) * pm; +% Make sure every gear is larger than the shaft that runs through it +check_planet_arr = [... + front_planet_s1_diameter, ... + front_planet_s2_diameter, ... + rear_planet_s1_diameter, ... + rear_planet_s2_diameter]; -g(7) = max(0, (planet_d2 / (0.45 * ring_d2)) - 1) * pm; +check_sun_arr = [... + front_sun_diameter, ... + rear_sun_diameter]; -if sun_d1 < sun_extrusion_d - g(8) = pm; -end -%% Sum constriants -val = sum(g); +% Count occurances where the constriant is violated +val = val + PENALTY * sum(check_planet_arr < planet_shaft_diameter); + +val = val + PENALTY * sum(check_sun_arr < sun_shaft_diameter); +% Enforce upper bounds on x +val = val + PENALTY * sum(x > maxX); diff --git a/evalFitness.m b/evalFitness.m new file mode 100644 index 0000000..99f0519 --- /dev/null +++ b/evalFitness.m @@ -0,0 +1,9 @@ +function fitnessVal = evalFitness(x, bearing_diameter, ... + planet_shaft_diameter, sun_shaft_diameter, maxX, GR_front_target, ... + GR_rear_target) + +phi = fval(x, GR_front_target, GR_rear_target); +penalty = constraint(x, bearing_diameter, planet_shaft_diameter, ... + sun_shaft_diameter, maxX); + +fitnessVal = phi + penalty; \ No newline at end of file diff --git a/explainX.m b/explainX.m new file mode 100644 index 0000000..083d296 --- /dev/null +++ b/explainX.m @@ -0,0 +1,66 @@ +function [front_diameters, front_teeth, ... + rear_diameters, rear_teeth, dp, frontGR, rearGR] = explainX(x) + +%% This function explains the design represented by x + +% X is formatted +% [ +% front_ring_teeth / 3 +% front_planet_s2_teeth +% front_sun_teeth / 3 +% rear_planet_s2_teeth +% rear_sun_teeth / 3 +% diametral_pitch +% ] + +%% Interpret x +front_ring_teeth = 3 * x(1); +front_planet_s2_teeth = x(2); +front_sun_teeth = 3 * x(3); +rear_planet_s2_teeth = x(4); +rear_sun_teeth = 3 * x(5); +dp = 12 + 4 * x(6); + + +%% Calculate derived variables +rear_ring_teeth = front_ring_teeth; + +front_planet_s1_teeth = front_ring_teeth - ... + front_planet_s2_teeth - front_sun_teeth; + +rear_planet_s1_teeth = rear_ring_teeth - ... + rear_planet_s2_teeth - rear_sun_teeth; + + +%% Calculate diameters +front_ring_diameter = front_ring_teeth / dp; +front_planet_s1_diameter = front_planet_s1_teeth / dp; +front_planet_s2_diameter = front_planet_s2_teeth / dp; +front_sun_diameter = front_sun_teeth / dp; + +rear_ring_diameter = rear_ring_teeth / dp; +rear_planet_s1_diameter = rear_planet_s1_teeth / dp; +rear_planet_s2_diameter = rear_planet_s2_teeth / dp; +rear_sun_diameter = rear_sun_teeth / dp; + + +%% Prepare outputs +front_diameters = [front_planet_s1_diameter, front_planet_s2_diameter, ... + front_sun_diameter, front_ring_diameter]; + +rear_diameters = [rear_planet_s1_diameter, rear_planet_s2_diameter, ... + rear_sun_diameter, rear_ring_diameter]; + +front_teeth = [front_planet_s1_teeth, front_planet_s2_teeth, ... + front_sun_teeth, front_ring_teeth]; + +rear_teeth = [rear_planet_s1_teeth, rear_planet_s2_teeth, ... + rear_sun_teeth, rear_ring_teeth]; + +frontGR = (front_planet_s1_teeth * front_ring_teeth) / ... + (front_sun_teeth * front_planet_s2_teeth) + 1; + +rearGR = (rear_planet_s1_teeth * rear_ring_teeth) / ... + (rear_sun_teeth * rear_planet_s2_teeth) + 1; + + diff --git a/fval.m b/fval.m index 1f901af..379bef4 100644 --- a/fval.m +++ b/fval.m @@ -1,13 +1,45 @@ -function val = fval(x, GR) -%% Planetary Gear Ratio Target Function -% Elliot Stockwell 1/12/21 + +function val = fval(x, GR_front_target, GR_rear_target) +%% Gear Ratio Objective Function % -% x : array of input variables (planet_d1, planet_d2, sun_d1, dp) -% GR : target gear ratio +% X is formatted +% [ +% front_ring_teeth / 3 +% front_planet_s2_teeth +% front_sun_teeth / 3 +% rear_planet_s2_teeth +% rear_sun_teeth / 3 +% diametral_pitch +% ] + + +%% Interpret x +front_ring_teeth = 3 * x(1); +front_planet_s2_teeth = x(2); +front_sun_teeth = 3 * x(3); +rear_planet_s2_teeth = x(4); +rear_sun_teeth = 3 * x(5); +dp = 12 + 4 * x(6); + + +%% Calculate derived variables +rear_ring_teeth = front_ring_teeth; + +front_planet_s1_teeth = front_ring_teeth - ... + front_planet_s2_teeth - front_sun_teeth; + +rear_planet_s1_teeth = rear_ring_teeth - ... + rear_planet_s2_teeth - rear_sun_teeth; + +%% Calculate gear ratios +front_gr = (front_planet_s1_teeth * front_ring_teeth) / ... + (front_sun_teeth * front_planet_s2_teeth) + 1; + +rear_gr = (rear_planet_s1_teeth * rear_ring_teeth) / ... + (rear_sun_teeth * rear_planet_s2_teeth) + 1; + -planet_d1 = x(1); -planet_d2 = x(2); -sun_d1 = x(3); -ring_d2 = planet_d1 + planet_d2 + sun_d1; +%% Calculate objective +val = (front_gr - GR_front_target) ^ 2; +val = val + (rear_gr - GR_rear_target) ^ 2; -val = (((planet_d1 * ring_d2 / sun_d1 / planet_d2) + 1) - GR)^2; \ No newline at end of file diff --git a/runAllOptimizations.m b/runAllOptimizations.m new file mode 100644 index 0000000..3f741c1 --- /dev/null +++ b/runAllOptimizations.m @@ -0,0 +1,120 @@ +%% Master running file for GA methods + +%% Define constants +bearing_diameter = 4.50; % in +sun_shaft_diameter = 1.23; % in +planet_shaft_diameter = 0.85; % in + +GR_front_target = 6.5; % dimensionless +GR_rear_target = 8.00; % dimensionless + +numGenerations = 200; +num_runs = 500; + +%% Define bounds +% +% X is formatted +% [ +% front_ring_teeth / 3 +% front_planet_s2_teeth +% front_sun_teeth / 3 +% rear_planet_s2_teeth +% rear_sun_teeth / 3 +% diametral_pitch +% ] + +% Remember some of these are / 3 +min_x = [... + 14, ... + 14, ... + 5, ... + 14, ... + 5, ... + 1]; + +max_x = [... + 35, ... + 48, ... + 35, ... + 48, ... + 35, ... + 4]; + + +%% Perform last calculations +vlb = min_x; + +vdiff = max_x - min_x; +bits = ceil(log(vdiff) / log(2)); + +vub = vlb + 2 .^ bits - 1; +maxX = max_x; + + +%% Run GA +data = zeros(num_runs, 20); + +fprintf('\nBeginning GA\n'); + +parfor i=1:num_runs + [FOPT, X] = runOneOptimization(vlb, vub, bits, bearing_diameter, ... + planet_shaft_diameter, sun_shaft_diameter, maxX, ... + GR_front_target, GR_rear_target, numGenerations) + + [front_diameters, front_teeth, rear_diameters, rear_teeth, ... + dp, frontGR, rearGR] = explainX(X); + + data(i, :) = [FOPT, front_diameters, front_teeth, ... + rear_diameters, rear_teeth, dp, frontGR, rearGR]; +end + + +%% Save Info +filename = 'ga_results.csv'; + +header = {'Badness', ... + 'Front Planet S1 Diameter (in)', 'Front Planet S2 Diameter (in)', ... + 'Front Sun Diameter (in)', 'Front Ring Diameter (in)', ... + 'Front Planet S1 Teeth', 'Front Planet S2 Teeth', ... + 'Front Sun Teeth', 'Front Ring Teeth', ... + 'Rear Planet S1 Diameter (in)', 'Rear Planet S2 Diameter (in)', ... + 'Rear Sun Diameter (in)', 'Rear Ring Diameter (in)', ... + 'Rear Planet S1 Teeth', 'Rear Planet S2 Teeth', ... + 'Rear Sun Teeth', 'Rear Ring Teeth', ... + 'Diametral Pitch (1/in)', 'Front GR', 'Rear GR'}; + +commaHeader = [header;repmat({','},1,numel(header))]; +commaHeader = commaHeader(:)'; +textHeader = cell2mat(commaHeader); + +fid = fopen(filename, 'w'); +fprintf(fid, '%s\n', textHeader(1:end-1)); +fclose(fid); + +data_sorted = sortrows(data); +dlmwrite(filename, data_sorted, '-append'); + + +%% Clean up +fprintf('\nDone.\n'); + +%% Create histogram + +feasible_region_front = [6.25 6.25 6.75 6.75 6.25]; +feasible_region_rear = [7.875 8.5 8.5 7.875 7.875]; + +figure(1); +hold on;grid on;grid minor; +colormap jet; +set(gcf, 'Position', [100 100 900 700]); +sct = scatter(data(:, end-1), data(:, end), 25, data(:, 1)); +trg = plot(GR_front_target, GR_rear_target); +region = plot(feasible_region_front, feasible_region_rear); +set(trg, 'Marker', '.', 'MarkerSize', 25, 'Color', 'r', ... + 'LineStyle', 'None'); +set(region, 'Color', 'k', 'LineStyle', '-.', 'LineWidth', 1); +title('Distribution of Results'); +xlabel('Front Gear Ratio');ylabel('Rear Gear Ratio'); +legend([sct, trg], 'GA Results', 'Target'); + + diff --git a/runOneOptimization.m b/runOneOptimization.m new file mode 100644 index 0000000..7645ebf --- /dev/null +++ b/runOneOptimization.m @@ -0,0 +1,20 @@ +function [FOPT, X] = runOneOptimization(vlb, vub, bits, bearing_diameter, ... + planet_shaft_diameter, sun_shaft_diameter, maxX, GR_front_target, ... + GR_rear_target, numGenerations) + +% Chromosome length +l = sum(bits); + +%% Define GA options +parain = zeros(14, 1); +parain(1) = 0; +parain(11) = 4 * l; +parain(13) = (l + 1) / (2 * l * parain(11)); +parain(14) = numGenerations; + +options = goptions(parain); + +[X, FOPT, ~, ~, ~, ~, ~] = GA550('evalFitness', ... + [], options, vlb, vub, bits, bearing_diameter, ... + planet_shaft_diameter, sun_shaft_diameter, maxX, GR_front_target, ... + GR_rear_target); \ No newline at end of file