diff --git a/generate_data.m b/generate_data.m index ab29ac4..ea76e90 100644 --- a/generate_data.m +++ b/generate_data.m @@ -465,18 +465,16 @@ display(sprintf(repmat('*', 1, 76))) display('Warning') display(sprintf(repmat('*', 1, 76))) - s = ['The closed loop system is not stable with these selected ', ... - 'controller gains. The Handling Quality Metric is not valid\n', ... - 'because the system is closed loop unstable. The simulations,', ... - 'will reflect the instabilities, but may be useful if\n', ... - 'marginally stable. It is best to provide improved gain guesses', ... - 'or supply the gains manually.\n']; + s = ['The closed loop system is not stable with these selected controller gains.\n', ... + 'The Handling Quality Metric is not valid because the system is closed loop\n', ... + 'unstable. The simulations will reflect the instabilities, but may be useful\n', ... + 'if marginally stable. It is best to provide improved gain guesses or supply\n', ... + 'the gains manually.']; display(sprintf(s)) display(sprintf(repmat('*', 1, 76))) display('Closed loop system eigenvalues:') eig(A) - display('Gains selected:') - display(settings.gains) + settings.gains else % if not unstable % only save gains if not user supplied and the neuro frequency is % default diff --git a/objective.m b/objective.m index fb101f8..86928ae 100644 --- a/objective.m +++ b/objective.m @@ -1,25 +1,66 @@ -function [peak_hqm] = objective(unknowns) +function [peak_hqm] = objective(unknowns, bicycle, speed) +% OBJECTIVE - Returns the maximum value of the HQM transfer function +% magnitude if the system is closed loop stable or an order of magnitude +% larger value based on the largest closed loop poles are if the system is +% closed loop unstable. +% +% Syntax: peak_hqm = objective(unknowns, bicycle, speed) +% +% Inputs: +% unknowns - Vector of the optimization parameters (c trail, w wheelbase, +% lam steer axis tilt, IFyy front wheel rotational inertia), +% size 4x1. +% bicycle - A char for the bicycle name, e.g. 'Benchmark', 'Pista', etc. +% speed - Scalar value for the design speed. +% +% Outputs: +% peak_hqm - Maximum scalar value of the HQM magnitude. -display('==========================================') - -unknowns +display(sprintf(repmat('=', 1, 79))) +display('Values passed into the objective function:') +display(sprintf('c: %1.5f ', unknowns(1))) +display(sprintf('w: %1.5f ', unknowns(2))) +display(sprintf('lam: %1.5f ', unknowns(3))) +display(sprintf('IFyy: %1.5f ', unknowns(4))) freqs = linspace(0.01, 40, 200); -par = par_text_to_struct('parameters/BenchmarkPar.txt'); +par = par_text_to_struct(['parameters/' bicycle 'Par.txt']); +% NOTE : Here to try for fun to see if different bicycle designs are needed +% for riding in low gravity. +% TODO : These should be moved to "optimal_bicycle.m". %par.g = 1.6; % acceleration due to gravity on the moon %par.g = 3.71; % acceleration due to gravity on mars -par.c = unknowns(1); -par.w = unknowns(2); -par.lam = unknowns(3); -par.rF = unknowns(4); -speed = 5.0; +par.c = unknowns(1); % trail +par.w = unknowns(2); % wheelbase +par.lam = unknowns(3); % steer axis tilt + +% Given an existing wheel with a rotational inertia, mass, and radius that +% approximately adhere to I=m*r^2, make sure that this relationship holds +% for the optimal front wheel inertias and that we only have to add mass to +% the existing wheel while keeping the radius constant or reduce the radius +% while keeping mass constant. These should be physically realizable barring +% that the radius doesn't get too tiny. +IFyy_opt = unknowns(4); + +if IFyy_opt < par.IFyy + % keep the mass the same but reduce the wheel radius + par.rF = sqrt(IFyy_opt ./ par.mF); +elseif IFyy_opt > par.IFyy + % keep radius of the wheel the same but add mass + par.mF = IFyy_opt ./ par.rF.^2; +end + +display(sprintf('rF: %1.5f ', par.rF)) +display(sprintf('mF: %1.5f ', par.mF)) + +par.IFyy = IFyy_opt; % front wheel rotational inertia [A, B, C, D] = whipple_pull_force_abcd(par, speed); -data = generate_data('Benchmark', speed, ... +data = generate_data(bicycle, speed, ... 'simulate', false, ... 'forceTransfer', {}, ... 'fullSystem', false, ... @@ -30,6 +71,11 @@ lateral_dev_loop = minreal(tf(data.closedLoops.Y.num, data.closedLoops.Y.den)); if ~isstable(lateral_dev_loop) + % TODO : It may be best to return NAN here as per the CMAES statement + % "An easy way to implement a hard non-linear constraint is to return + % NaN. Then, this function evaluation is not counted and a newly sampled + % point is tried immediately." Right now it returns a large number + % relative to the target low HQM values. peak_hqm = max(10, 100 * max(real(pole(lateral_dev_loop)))); else num = data.handlingMetric.num; @@ -38,4 +84,6 @@ peak_hqm = max(mag); end +display('Value of the objective function:') peak_hqm +display(sprintf(repmat('=', 1, 79))) diff --git a/optimal_bicycle.m b/optimal_bicycle.m index 650a01e..1260b70 100644 --- a/optimal_bicycle.m +++ b/optimal_bicycle.m @@ -1,18 +1,38 @@ -par = par_text_to_struct('parameters/BenchmarkPar.txt'); +% Searches for a set of four parameters (c, w, lam, IFyy) that minimizes the +% peak HQM magnitude for a specific design speed. + +% Set the design speed and the base bicycle here: +speed = 3.0; % m/s +bicycle = 'Pista'; + +par = par_text_to_struct(['parameters/' bicycle 'Par.txt']); guess = zeros(4, 1); guess(1) = par.c; guess(2) = par.w; guess(3) = par.lam; -guess(4) = par.rF; +guess(4) = par.IFyy; -% wheelbase should always accomdate the mass center +% wheelbase should always accomdate the mass center (bicycle doesn't pitch +% foward) min_wheelbase = (par.mH * par.xH + par.mB * par.xB) / (par.mH + par.mB); -opts.LBounds = [-inf; min_wheelbase; -pi/2; 1e-10]; -opts.UBounds = [inf; inf; pi/2; inf]; +fprintf('The minimum possible wheelbase is %1.4f\n', min_wheelbase); + +opts.LBounds = [-inf; % c + min_wheelbase; % w + -pi/2; % lam + 3e-5]; % IFyy + +opts.UBounds = [inf; % c + inf; % w + pi/2; % lam + inf]; % IFyy -sigma = [0.5; 3.0; 0.3 * pi; 0.2]; +sigma = [0.5; + 3.0; + 0.3 * pi; + 0.07]; %sigma = sqrt(var(guess')'); -[optimal_par, hqm] = cmaes('objective', guess, sigma, opts); +[optimal_par, hqm] = cmaes('objective', guess, sigma, opts, bicycle, speed);