Files
linguist/samples/Matlab/create_ieee_paper_plots.m
Jason Moore 9bb230d7c8 Added matlab code samples.
All of these code samples currently are mis-identified in my repositories. I'm
donating them to the cause.
2013-01-30 13:12:45 -08:00

1120 lines
33 KiB
Matlab

function create_ieee_paper_plots(data, rollData)
% Creates all of the figures for the IEEE paper.
%
% Parameters
% ----------
% data : structure
% A structure contating the data from generate_data.m for all of the bicycles
% and speeds for the IEEE paper.
% rollData : structure
% The data for a single bicycle at a single speed with roll torque as the
% input.
global goldenRatio
% used for figure width to height ratio
goldenRatio = (1 + sqrt(5)) / 2;
% create a plot directory if one doesn't already exist
if exist('plots/', 'dir') ~= 7
mkdir('plots/')
end
% Define some linestyles and colors for each of the six bicycles
linestyles = {'-', '-', '-.', ...
'--', '-.', '--'};
colors = {'k', ...
[0.5, 0.5, 0.5], ...
[0.5, 0.5, 0.5], ...
'k', ...
'k', ...
[0.5, 0.5, 0.5]};
loop_shape_example(data.Benchmark.Medium, 'Steer')
loop_shape_example(rollData, 'Roll')
plot_io_roll(rollData, 'Distance')
plot_io_roll(rollData, 'Time')
open_loop_all_bikes(data, linestyles, colors)
handling_all_bikes(data, rollData, linestyles, colors)
path_plots(data, linestyles, colors)
var = {'delta', 'phi', 'psi', 'Tdelta'};
io = {'output', 'output', 'output', 'input'};
typ = {'Distance', 'Time'};
for i = 1:length(var)
for j = 1:length(typ)
plot_io(var{i}, io{i}, typ{j}, data, linestyles, colors)
end
end
phase_portraits(data.Benchmark.Medium)
eigenvalues(data, linestyles, colors)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function loop_shape_example(bikeData, input)
% Creates the example loop shaping for the bicycle at medium speed.
%
% Parameters
% ----------
% bikeData : structure
% Contains data for a single bicycle at a single speed.
% input : string
% 'Steer' or 'Roll' depending on what input was used to control the bicycle.
global goldenRatio
% closed loop bode plots
figure()
figWidth = 5.0;
figHeight = figWidth / goldenRatio;
set(gcf, ...
'Color', [1, 1, 1], ...
'PaperOrientation', 'portrait', ...
'PaperUnits', 'inches', ...
'PaperPositionMode', 'manual', ...
'OuterPosition', [424, 305 - 50, 518, 465], ...
'PaperPosition', [0, 0, figWidth, figHeight], ...
'PaperSize', [figWidth, figHeight])
freq = {0.1, 20.0};
hold all
closedLoops = bikeData.closedLoops;
% make sure all bode plots display 'rad/s' instead of 'rad/sec'
bops = bodeoptions;
bops.FreqUnits = 'rad/s';
if strcmp(input, 'Steer')
linestyles = {'', '', '-.', '-.', '-', '-.', '-.', '-'};
gray = [0.6, 0.6, 0.6];
colors = {'k', 'k', 'k', gray, 'k', 'k', gray, 'k'};
% the closed delta loop
deltaNum = closedLoops.Delta.num;
deltaDen = closedLoops.Delta.den;
bodeplot(tf(deltaNum, deltaDen), freq, bops);
% a typical neuromuscular model
neuroNum = 2722.5;
neuroDen = [1, 13.96, 311.85, 2722.5];
bodeplot(tf(neuroNum, neuroDen), freq, bops);
whichLines = 5:-1:3;
elseif strcmp(input, 'Roll')
linestyles = {'', '', '-', '-'};
colors = {'k', 'k', 'k', 'k'};
whichLines = 4:-1:2;
else
error('Bad input, use Steer or Roll')
end
% the closed phi dot loop
phiDotNum = closedLoops.PhiDot.num;
phiDotDen = closedLoops.PhiDot.den;
closedBode = bodeplot(tf(phiDotNum, phiDotDen), freq, bops);
hold off
% clean it up
opts = getoptions(closedBode);
if strcmp(input, 'Steer')
opts.YLim = {[-45, 20], [-360, 90]};
else
opts.YLim = {[-30, 10], [-180, 90]};
end
opts.PhaseMatching = 'on';
opts.PhaseMatchingValue = 0;
opts.Title.String = '';
setoptions(closedBode, opts)
% find all the lines in the current figure
lines = findobj(gcf, 'type', 'line');
for i = 3:length(lines)
set(lines(i), 'LineStyle', linestyles{i}, ...
'Color', colors{i}, ...
'LineWidth', 2.0)
end
% there seems to be a bug such that the xlabel is too low, this is a hack to
% get it to work
raise = 0.05;
plotAxes = findobj(gcf, 'type', 'axes');
set(plotAxes, 'XColor', 'k', 'YColor', 'k')
curPos1 = get(plotAxes(1), 'Position');
curPos2 = get(plotAxes(2), 'Position');
set(plotAxes(1), 'Position', curPos1 + [0, raise, 0, 0])
set(plotAxes(2), 'Position', curPos2 + [0, raise, 0, 0])
xLab = get(plotAxes(1), 'Xlabel');
set(xLab, 'Units', 'normalized')
set(xLab, 'Position', get(xLab, 'Position') + [0, raise + 0.05, 0])
% make the tick labels smaller
set(plotAxes, 'Fontsize', 8)
if strcmp(input, 'Steer')
legWords = {'$\delta$ Loop',
'Neuromuscular model from [27]',
'$\dot{\phi}$ Loop'};
elseif strcmp(input, 'Roll')
legWords = {'$\dot{\phi}$ Loop'};
end
closeLeg = legend(lines(whichLines), ...
legWords, ...
'Location', 'Southwest', ...
'Interpreter', 'Latex', ...
'Fontsize', 8);
% add the annotation showing a 10 dB peak
if strcmp(input, 'Steer')
axes(plotAxes(2))
db1 = text(2.7, 5.0, '~10dB');
db2 = text(2.5, -10.0, '~10dB');
set([db1, db2], 'Fontsize', 8)
dArrow1 = annotation('doublearrow', ...
[0.7, 0.7], ...
[0.755 + raise, 0.818 + raise]);
annotation('line', [0.69, 0.87], [0.818 + raise, 0.818 + raise])
dArrow2 = annotation('doublearrow', ...
[0.685, 0.685], ...
[0.665 + raise, 0.725 + raise]);
annotation('line', [0.675, 0.87], [0.725 + raise, 0.725 + raise])
set([dArrow1, dArrow2], 'Head1width', 3, 'Head1length', 3, ...
'Head2width', 3, 'Head2length', 3)
else
axes(plotAxes(2))
db1 = text(0.67, -3.7, '~10dB');
set(db1, 'Fontsize', 8)
dArrow = annotation('doublearrow', ...
[0.5, 0.5], ...
[0.697 + raise, 0.795 + raise]);
set(dArrow, 'Head1width', 3, 'Head1length', 3, ...
'Head2width', 3, 'Head2length', 3)
annotation('line', [0.49, 0.75], [0.795 + raise, 0.795 + raise])
end
filename = ['benchmark' input 'Closed'];
pathToFile = ['plots' filesep filename];
print(gcf, '-deps2c', '-loose', [pathToFile '.eps'])
fix_ps_linestyle([pathToFile '.eps'])
% open loop plots
figure()
set(gcf, ...
'Color', [1, 1, 1], ...
'PaperOrientation', 'portrait', ...
'PaperUnits', 'inches', ...
'PaperPositionMode', 'manual', ...
'PaperPosition', [0, 0, figWidth, figHeight], ...
'PaperSize', [figWidth, figHeight])
openLoops = bikeData.openLoops;
hold all
num = openLoops.Phi.num;
den = openLoops.Phi.den;
bodeplot(tf(num, den), freq, bops);
num = openLoops.Psi.num;
den = openLoops.Psi.den;
bodeplot(tf(num, den), freq, bops);
num = openLoops.Y.num;
den = openLoops.Y.den;
openBode = bodeplot(tf(num, den), freq, bops);
hold off
% clean it up
opts = getoptions(openBode);
opts.Title.String = '';
opts.YLim = {[-80, 20], [-540, -60]};
opts.PhaseMatching = 'on';
opts.PhaseMatchingValue = 0;
setoptions(openBode, opts)
% find all the lines in the current figure
lines = findobj(gcf, 'type', 'line');
linestyles = {'', '', '-.', '-', '--', '-.', '-', '--'};
for i = 3:length(lines)
set(lines(i), 'LineStyle', linestyles{i}, ...
'Color', 'k', ...
'LineWidth', 2.0)
end
plotAxes = findobj(gcf, 'type', 'axes');
set(plotAxes, 'Fontsize', 8, 'XColor', 'k', 'YColor', 'k')
closeLeg = legend(lines(8:-1:6), ...
{'$\phi$ Loop', '$\psi$ Loop','$y$ Loop'}, ...
'Location', 'Southwest', ...
'Interpreter', 'Latex');
% add zero crossing lines
%axes(plotAxes(1))
%line([0.1, 20], [-180, -180], 'Color', 'k')
axes(plotAxes(2))
line([0.1, 20], [0, 0], 'Color', 'k')
% add some lines and labels for the cross over frequencies
if strcmp(input, 'Steer')
wc = 2;
wShift = [0.42, 0.35, 0.175];
else strcmp(input, 'Roll')
wc = 1.5;
wShift = [0.31, 0.26, 0.1325];
end
axes(plotAxes(2))
hold on
gray = [0.5, 0.5, 0.5];
line([wc, wc], [-40, 0], 'Color', gray)
text(wc - wShift(1), -43, ['$\omega_c=' num2str(wc) '$'], ...
'Interpreter', 'Latex', 'Fontsize', 8)
line([wc / 2, wc / 2], [-30, 0], 'Color', gray)
text(wc / 2 - wShift(2), -33, ['$\omega_c/2=' num2str(wc / 2) '$'], ...
'Interpreter', 'Latex', 'Fontsize', 8)
line([wc / 4, wc / 4], [-20, 0], 'Color', gray)
text(wc / 4 - wShift(3), -23, ['$\omega_c/4=' num2str(wc / 4) '$'], ...
'Interpreter', 'Latex', 'Fontsize', 8)
hold off
curPos1 = get(plotAxes(1), 'Position');
curPos2 = get(plotAxes(2), 'Position');
set(plotAxes(1), 'Position', curPos1 + [0, raise, 0, 0])
set(plotAxes(2), 'Position', curPos2 + [0, raise, 0, 0])
xLab = get(plotAxes(1), 'Xlabel');
set(xLab, 'Units', 'normalized')
set(xLab, 'Position', get(xLab, 'Position') + [0, raise + 0.05, 0])
filename = ['benchmark' input 'Open.eps'];
pathToFile = ['plots' filesep filename];
print(gcf, '-deps2c', '-loose', pathToFile)
fix_ps_linestyle(pathToFile)
% handling qualities plot
num = bikeData.handlingMetric.num;
den = bikeData.handlingMetric.den;
w = linspace(0.01, 20, 200);
[mag, phase, freq] = bode(tf(num, den), w);
figure()
set(gcf, ...
'Color', [1, 1, 1], ...
'PaperOrientation', 'portrait', ...
'PaperUnits', 'inches', ...
'PaperPositionMode', 'manual', ...
'PaperPosition', [0, 0, figWidth, figHeight], ...
'PaperSize', [figWidth, figHeight])
hold on
metricLine = plot(freq, mag(:)', 'k-', 'Linewidth', 2.0);
level1 = line([0, 20], [5, 5]);
level2 = line([0, 20], [8, 8]);
set(level1, 'Color', 'k', 'Linestyle', '--', 'Linewidth', 2.0)
set(level2, 'Color', 'k', 'Linestyle', '--', 'Linewidth', 2.0)
ylim([0, 10]);
ylabel('Handling Quality Metric')
xlabel('Frequency (rad/s)')
text(3, 3, 'Level 1')
text(3, 6.5, 'Level 2')
text(3, 9, 'Level 3')
box on
filename = ['benchmark' input 'Handling.eps'];
pathToFile = ['plots' filesep filename];
print(gcf, '-deps2', '-loose', pathToFile)
fix_ps_linestyle(pathToFile)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function open_loop_all_bikes(data, linestyles, colors)
% Creates open loop Bode plots of all the bikes.
global goldenRatio
bikes = fieldnames(data);
figure()
figWidth = 5.0;
figHeight = figWidth / goldenRatio;
set(gcf, ...
'Color', [1, 1, 1], ...
'PaperOrientation', 'portrait', ...
'PaperUnits', 'inches', ...
'PaperPositionMode', 'manual', ...
'PaperPosition', [0, 0, figWidth, figHeight], ...
'PaperSize', [figWidth, figHeight])
freq = {0.1, 20.0};
% make sure all bode plots display 'rad/s' instead of 'rad/sec'
bops = bodeoptions;
bops.FreqUnits = 'rad/s';
hold all
for i = 2:length(bikes)
num = data.(bikes{i}).Medium.openLoops.Phi.num;
den = data.(bikes{i}).Medium.openLoops.Phi.den;
openBode = bodeplot(tf(num, den), freq, bops);
end
hold off
% clean it up
opts = getoptions(openBode);
%opts.Title.String = '$\phi$ Open Loop Bode Diagrams at 5 m/s';
opts.Title.String = '';
%opts.Title.Interpreter = 'Latex';
opts.YLim = {[-30, 10], [-540, -90]};
opts.PhaseMatching = 'on';
opts.PhaseMatchingValue = 0;
setoptions(openBode, opts)
% find all the lines in the current figure
plotAxes = findobj(gcf, 'type', 'axes');
magLines = findobj(plotAxes(2), 'type', 'line');
phaseLines = findobj(plotAxes(1), 'type', 'line');
for i = 2:length(magLines)
set(magLines(i), ...
'LineStyle', linestyles{i - 1}, ...
'Color', colors{i - 1}, ...
'LineWidth', 1.0)
set(phaseLines(i), ...
'LineStyle', linestyles{i - 1}, ...
'Color', colors{i - 1}, ...
'LineWidth', 1.0)
end
closeLeg = legend(magLines(2:7), ...
{'1', '2', '3', '4', '5', '6'}, ...
'Location', 'Southwest', ...
'Fontsize', 8);
set(plotAxes, 'YColor', 'k', 'XColor', 'k', 'Fontsize', 8)
% add a zero lines
axes(plotAxes(1))
line([0.1, 20], [-180, -180], 'Color', 'k')
axes(plotAxes(2))
line([0.1, 20], [0, 0], 'Color', 'k')
% raise the axes cause the xlabel is cut off
raise = 0.05;
curPos1 = get(plotAxes(1), 'Position');
curPos2 = get(plotAxes(2), 'Position');
set(plotAxes(1), 'Position', curPos1 + [0, raise, 0, 0])
set(plotAxes(2), 'Position', curPos2 + [0, raise, 0, 0])
xLab = get(plotAxes(1), 'Xlabel');
set(xLab, 'Units', 'normalized')
set(xLab, 'Position', get(xLab, 'Position') + [0, raise + 0.05, 0])
filename = 'openBode.eps';
pathToFile = ['plots' filesep filename];
print(pathToFile, '-deps2c', '-loose')
fix_ps_linestyle(pathToFile)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function handling_all_bikes(data, rollData, linestyles, colors)
% Creates handling quality metric for all bikes.
%
% Parameters
% ----------
% data : structure
% Contains data for all bikes a the three speeds for steer input.
% rollData : structure
% Contains the data for the benchmark bike with roll input at medium speed.
% linestyles : cell array
% Linestyle strings, one for each of the six bikes.
% colors : cell array
% Colorspecs for each of the six bikes.
global goldenRatio
bikes = fieldnames(data);
figure()
figWidth = 5.0;
figHeight = figWidth / goldenRatio;
set(gcf, ...
'Color', [1, 1, 1], ...
'PaperOrientation', 'portrait', ...
'PaperUnits', 'inches', ...
'PaperPositionMode', 'manual', ...
'PaperPosition', [0, 0, figWidth, figHeight], ...
'PaperSize', [figWidth, figHeight])
w = linspace(0.01, 20, 200);
speedNames = fieldnames(data.Browser);
fillColors = {[0.82, 0.82, 0.82]
[0.68, 0.68, 0.68]
[0.95, 0.95, 0.95]};
hold all
% plot the background area for each family of curves
for j = 1:length(speedNames)
% get the max values for the set of curves
magnitudes = zeros(length(w), length(bikes) - 1);
for i = 2:length(bikes)
num = data.(bikes{i}).(speedNames{j}).handlingMetric.num;
den = data.(bikes{i}).(speedNames{j}).handlingMetric.den;
[mag, phase, freq] = bode(tf(num, den), w);
magnitudes(:, i - 1) = mag(:)';
end
maxMag = max(magnitudes, [], 2);
% fill the area under the curve
area(freq, maxMag, ...
'Facecolor', fillColors{j}, ...
'Edgecolor', 'none')
end
% this makes sure that the edges of the fill area don't cover the axes
set(gca, 'Layer', 'top')
% plot the actual curves
for j = 1:length(speedNames)
metricLines = zeros(length(bikes) - 1, 1);
for i = 2:length(bikes)
num = data.(bikes{i}).(speedNames{j}).handlingMetric.num;
den = data.(bikes{i}).(speedNames{j}).handlingMetric.den;
[mag, phase, freq] = bode(tf(num, den), w);
metricLines(i - 1) = plot(freq, mag(:)', ...
'Color', colors{i - 1}, ...
'Linestyle', linestyles{i - 1}, ...
'Linewidth', 2.0);
end
end
% add the roll input bike
num = rollData.handlingMetric.num;
den = rollData.handlingMetric.den;
[mag, phase, freq] = bode(tf(num, den), w);
rollLine = plot(freq, mag(:)', 'k', 'Linewidth', 2.0, 'Linestyle', ':');
hold off
% move the roll input line down so it shows on the legend
chil = get(gca, 'Children');
legLines = [chil(end:-1:14)', rollLine];
legend(legLines, [{'2.5 m/s', '5.0 m/s', '7.5 m/s'}, ...
{'1', '2', '3', '4', '5', '6', 'Hands-free @ 5 m/s'}], ...
'Fontsize', 8)
ylim([0, 20]);
level1 = line([0, 20], [5, 5]);
level2 = line([0, 20], [8, 8]);
set(level1, 'Color', 'k', 'Linestyle', '--', 'Linewidth', 1.0)
set(level2, 'Color', 'k', 'Linestyle', '--', 'Linewidth', 1.0)
ylabel('Handling Quality Metric')
xlabel('Frequency (rad/s)')
text(3.1, 4.3, 'Level 1')
text(1.9, 6.5, 'Level 2')
text(3, 15, 'Level 3')
box on
set(gca, 'YColor', 'k')
filename = 'handling.eps';
pathToFile = ['plots' filesep filename];
print(pathToFile, '-deps2c', '-loose')
fix_ps_linestyle(pathToFile)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function path_plots(data, linestyles, colors)
% Creates a plot of the path tracking for all bikes at all speeds.
global goldenRatio
bikes = fieldnames(data);
speedNames = fieldnames(data.Browser);
figure()
figWidth = 5.0;
figHeight = figWidth / goldenRatio;
set(gcf, ...
'Color', [1, 1, 1], ...
'PaperOrientation', 'portrait', ...
'PaperUnits', 'inches', ...
'PaperPositionMode', 'manual', ...
'PaperPosition', [0, 0, figWidth, figHeight], ...
'PaperSize', [figWidth, figHeight])
hold all
% shifts the paths by this many meters
shift = [0, 15, 35];
for j = 1:length(speedNames)
time = data.(bikes{2}).(speedNames{j}).time;
path = data.(bikes{2}).(speedNames{j}).path;
speed = data.(bikes{2}).(speedNames{j}).speed;
plot(time * speed + shift(j), -path * j, 'k-')
for i = 2:length(bikes)
x = data.(bikes{i}).(speedNames{j}).outputs(:, 17);
x = x + shift(j);
y = data.(bikes{i}).(speedNames{j}).outputs(:, 18);
plot(x, -y * j, ...
'Linestyle', linestyles{i - 1}, ...
'Color', colors{i - 1}, ...
'Linewidth', 0.75)
end
[minPath, minPathI] = min(-path * j);
dis = time * speed + shift(j);
lab = sprintf('%1.1f m/s', speed);
text(dis(minPathI) - 15, minPath - 0.4, lab)
end
hold off
% change the y tick labels to positive and to reflect the 2 meter width
set(gca, 'YTick', [-7, -6, -4, -2, 0, 1])
set(gca, 'YTickLabel', {'', '2', '2', '2', '0', ''})
xlim([30 200])
box on
legend(['Path', {'1', '2', '3', '4', '5', '6'}], ...
'Fontsize', 8, 'Location', 'Southeast')
xlabel('Distance (m)')
ylabel('Lateral Deviation (m)')
filename = 'paths.eps';
pathToFile = ['plots' filesep filename];
print(pathToFile, '-deps2c', '-loose')
fix_ps_linestyle(pathToFile)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function plot_io(variable, io, xAxis, data, linestyles, colors)
% Creates a plot of the time histories of a particular output or input variable
% for three different speeds with either time or distance on the x axis.
%
% Parameters
% ----------
% variable : string
% The name of the variable you'd like to plot.
% io : string
% 'input' for input and 'output' for output.
% data : structure
% Data for a set of bicycles, the first being the benchmark bicycle.
% xAxis : string
% 'Distance' or 'Time' on the x axis.
% linestyles : cell array
% An array of linestyle types, one for each bicycle.
% colors : cell array
% An array of colors, one for each bicycle.
global goldenRatio
if strcmp(io, 'input')
names = {'Tphi',
'Tdelta',
'F'};
prettyNames = {'$T_\phi$',
'$T_\delta$',
'$F$'};
units = {'(N-m)',
'(N-m)',
'(N)'};
elseif strcmp(io, 'output')
names = {'xP',
'yP',
'psi',
'phi',
'thetaP',
'thetaR',
'delta',
'thetaF',
'xPDot',
'ypDot',
'psiDot',
'phiDot',
'thetaPDot',
'thetaRDot',
'deltaDot',
'thetaFDot',
'xQ',
'yQ'};
prettyNames = {'$x_P$',
'$y_P$',
'$\psi$',
'$\phi$',
'$\theta_P$',
'$\theta_R$',
'$\delta$',
'$\theta_F$',
'$\dot{x}_P$',
'$\dot{y}_P$',
'$\dot{\psi}$',
'$\dot{\phi}$',
'$\dot{\theta}_P$',
'$\dot{\theta}_R$',
'$\dot{\delta}$',
'$\dot{\theta}_F$',
'$x_Q$',
'$y_Q$'};
units = {'(m)',
'(m)',
'(rad)',
'(rad)',
'(rad)',
'(rad)',
'(rad)',
'(rad)',
'(m/s)',
'(m/s)',
'(rad/s)',
'(rad/s)',
'(rad/s)',
'(rad/s)',
'(rad/s)',
'(rad/s)',
'(m)',
'(m)'};
else
error('Please choose i or o')
end
index = find(ismember(names, variable) == 1);
bikes = fieldnames(data);
speedNames = fieldnames(data.Browser);
figure()
figWidth = 5.0;
figHeight = figWidth / goldenRatio;
set(gcf, ...
'Color', [1, 1, 1], ...
'PaperOrientation', 'portrait', ...
'PaperUnits', 'inches', ...
'PaperPositionMode', 'manual', ...
'PaperPosition', [0, 0, figWidth, figHeight], ...
'PaperSize', [figWidth, figHeight])
% find the maximum value of the variable
maxValue = 0;
for i = 2:length(bikes)
for j = 1:length(speedNames)
oneSpeed = data.(bikes{i}).(speedNames{j});
history = oneSpeed.([io 's'])(:, index);
if max(history) > maxValue
maxValue = max(history);
end
end
end
m = round(maxValue * 100) / 100;
pad = 0.15 * m;
yShift = [0, 2 * (m + pad), 4 * (m + pad)];
% shifts the paths by this many meters along the x axis
xShift = [0, 15, 35];
hold all
for j = 1:length(speedNames)
for i = 2:length(bikes)
oneSpeed = data.(bikes{i}).(speedNames{j});
time = oneSpeed.time;
speed = oneSpeed.speed;
distance = time * speed + xShift(j);
% time history of the value
history = oneSpeed.([io 's'])(:, index) + yShift(j);
if strcmp(xAxis, 'Distance')
xData = distance;
textX = 165;
elseif strcmp(xAxis, 'Time')
xData = time;
textX = 2;
else
error('Choose Time or Distance, no other')
end
plot(xData, history, ...
'Linestyle', linestyles{i - 1}, ...
'Color', colors{i - 1}, ...
'Linewidth', 0.75)
end
% add labels for the speeds
text(textX, yShift(j) + 4 * pad, [num2str(speed) ' m/s'])
end
ylim([-m - pad, yShift(3) + m + pad])
set(gca, 'YTick', ...
[-m, yShift(1), m, ...
yShift(2) - m, yShift(2), yShift(2) + m, ...
yShift(3) - m, yShift(3), yShift(3) + m])
ticks = {num2str(-m), '0', num2str(m)};
set(gca, 'YTickLabel', [ticks, ticks ticks])
if strcmp(xAxis, 'Distance')
xlabel('Distance (m)')
xLimits = [35, 190];
xlim(xLimits)
loc = 'Northwest';
else
xlabel('Time (s)')
xLimits = [0, 50];
xlim(xLimits)
loc = 'Northeast';
end
l1 = line(xLimits, [yShift(1) + m + pad, yShift(1) + m + pad]);
l2 = line(xLimits, [yShift(2) + m + pad, yShift(2) + m + pad]);
set([l1, l2], 'Color', 'k')
hold off
set(gca, 'Fontsize', 8)
first = [prettyNames{index} ' ' units{index}];
ylabel(first, 'Interpreter', 'Latex')
box on
legend({'1', '2', '3', '4', '5', '6'}, 'Fontsize', 8, 'Location', loc)
% if it is the steer angle plot for distance, add a magnifier for the
% countersteer
if strcmp(variable, 'delta') && strcmp(xAxis, 'Distance')
% Specify the position and the size of the rectangle
x_r = 37; y_r = 0; w_r = 4; h_r = 0.01;
rectangle('Position', [x_r-w_r/2, y_r-h_r/2, w_r, h_r], ...
'EdgeColor', 'k');
% Specify the position and the size of the 2. axis
x_a = 0.2; y_a = 0.29; w_a = 0.15; h_a = w_a * h_r / w_r * 20 / 0.05;
ax = axes('Units', 'Normalized', ...
'Position', [x_a, y_a, w_a, h_a], ...
'XTick', [], ...
'YTick', [], ...
'Box', 'on', ...
'LineWidth', 0.5, ...
'Color', 'w');
hold on
j = 1;
for i = 2:length(bikes)
oneSpeed = data.(bikes{i}).(speedNames{j});
time = oneSpeed.time;
speed = oneSpeed.speed;
distance = time * speed + xShift(j);
% time history of the value
history = oneSpeed.([io 's'])(:, index) + yShift(j);
plot(distance, history, ...
'Linestyle', linestyles{i - 1}, ...
'Color', colors{i - 1}, ...
'Linewidth', 0.75)
end
hold off
axis([x_r-w_r/2, x_r+w_r/2, y_r-h_r/2, y_r+h_r/2]);
text(35.3, -0.003, 'Countersteer', 'Fontsize', 8)
% bottom left
annotation('line', [x_a, 0.129], [y_a, 0.235])
% top left
annotation('line', [x_a, 0.132], [y_a + h_a, 0.26])
% bottom right
annotation('line', [x_a + w_a, 0.15], [y_a, 0.235])
% top right
annotation('line', [x_a, 0.15], [y_a + 0.02, 0.26])
end
% save the file
filename = [variable xAxis '.eps'];
print(['plots' filesep filename], '-deps2c', '-loose')
fix_ps_linestyle(['plots' filesep filename])
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function plot_io_roll(rollData, xAxis)
global goldenRatio
% closed loop bode plots
figure()
figWidth = 5.0;
figHeight = figWidth / goldenRatio;
set(gcf, ...
'Color', [1, 1, 1], ...
'PaperOrientation', 'portrait', ...
'PaperUnits', 'inches', ...
'PaperPositionMode', 'manual', ...
'PaperPosition', [0, 0, figWidth, figHeight], ...
'PaperSize', [figWidth, figHeight])
speed = rollData.speed;
time = rollData.time;
path = rollData.path;
frontWheel = rollData.outputs(:, 18);
rollAngle = rollData.outputs(:, 4);
steerAngle = rollData.outputs(:, 7);
rollTorque = rollData.inputs(:, 1);
% plot the path
subplot(2, 1, 1)
hold all
if strcmp(xAxis, 'Distance')
plot(speed * time, -path, 'k-', 'Linewidth', 1.0)
plot(speed * time, -frontWheel, 'k:', 'Linewidth', 1.0)
xlabel('Distance (m)')
xlim([30, 150])
elseif strcmp(xAxis, 'Time')
plot(time, -path, 'k-', 'Linewidth', 1.0)
plot(time, -frontWheel, 'k:', 'Linewidth', 1.0)
xlabel('Time (s)')
xlim([30 / speed, 150 / speed])
else
error('Bad xAxis, choose Distance or Time')
end
hold off
box on
ylabel('Lateral Deviation (m)')
ylim([-2.2, 0.2])
set(gca, 'YTickLabel', {'2', '1', '0'})
legend({'Path'}, ...
'Interpreter', 'Latex', ...
'Fontsize', 8, ...
'Location', 'Southeast')
subplot(2, 1, 2)
hold all
if strcmp(xAxis, 'Distance')
plot(speed * time, rollAngle, 'k-', 'Linewidth', 1.0)
[ax, h1, h2] = plotyy(speed * time, steerAngle, speed * time, rollTorque);
xlabel('Distance (m)')
xlim(ax(1), [30, 150])
xlim(ax(2), [30, 150])
elseif strcmp(xAxis, 'Time')
plot(time, rollAngle, 'k-', 'Linewidth', 1.0)
[ax, h1, h2] = plotyy(time, steerAngle, time, rollTorque);
xlabel('Time (s)')
xlim(ax(1), [30 / speed, 150 / speed])
xlim(ax(2), [30 / speed, 150 / speed])
else
error('Bad xAxis, choose Distance or Time')
end
hold off
box on
set(get(ax(1), 'Ylabel'), ...
'String', 'Angle (rad)', ...
'Color', 'k')
set(get(ax(2), 'Ylabel'), ...
'String', 'Torque (N-m)', ...
'Color', 'k')
set(ax, 'YColor', 'k', 'Fontsize', 8)
set(h1, 'Linestyle', '--', 'Color', 'k', 'Linewidth', 1.0)
set(h2, 'Linestyle', '-.', 'Color', 'k', 'Linewidth', 1.0)
legend({'$\phi$', '$\delta$', '$T_\phi$'}, ...
'Interpreter', 'Latex', ...
'Fontsize', 8, ...
'Location', 'Northeast')
if strcmp(xAxis, 'Distance')
axes(ax(2))
% magnifier rectangle
x_r = 38; y_r = 0; w_r = 10; h_r = 0.3;
rectangle('Position', [x_r-w_r/2, y_r-h_r/2, w_r, h_r], ...
'EdgeColor', 'k');
% magnify it
x_a = 0.14; y_a = 0.432; w_a = 0.28; h_a = 0.12;
inset = axes('Units', 'Normalized', ...
'Position', [x_a, y_a, w_a, h_a], ...
'Box', 'on', ...
'LineWidth', 0.5, ...
'Color', [0.8, 0.8, 0.8]);
hold on
plot(inset, time, rollAngle, 'k-', 'Linewidth', 1.0)
[ax, h1, h2] = plotyy(inset, time, steerAngle, time, rollTorque);
set(ax, 'XTick', [], ...
'YTick', [], ...
'YColor', 'k')
set(h1, 'Linestyle', '--', 'Color', 'k', 'Linewidth', 1.0)
set(h2, 'Linestyle', '-.', 'Color', 'k', 'Linewidth', 1.0)
axis(ax(1), [7, 8.5, -0.0025, 0.0025])
axis(ax(2), [7, 8.5, -0.0025 / 0.02, 0.0025 / 0.02])
% draw some lines connecting the corners
annotation('line', [x_a, 0.149], [y_a, 0.287])
annotation('line', [x_a + w_a, 0.213], [y_a, 0.287])
annotation('textbox', [0.26, 0.47, 0.1, 0.02], ...
'String', 'Countersteer', ...
'Fontsize', 8, ...
'Edgecolor', 'none')
end
filename = ['roll' xAxis '.eps'];
print(gcf, ['plots' filesep filename], '-deps2c', '-loose')
fix_ps_linestyle(['plots' filesep filename])
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function phase_portraits(bikeData)
% Creates four phase portrait plots to demonstrate the effects on the phase
% portraits when adjusting the gains.
global goldenRatio
figure()
figWidth = 6.0;
figHeight = figWidth / goldenRatio;
set(gcf, ...
'Color', [1, 1, 1], ...
'PaperOrientation', 'portrait', ...
'PaperUnits', 'inches', ...
'PaperPositionMode', 'manual', ...
'OuterPosition', [424, 305 - 50, 518, 465], ...
'PaperPosition', [0, 0, figWidth, figHeight], ...
'PaperSize', [figWidth, figHeight])
% this is the gain multiplier for the non-nominal plot
gainChanges = [1.2, 1, 1, 1, 1;
1, 1.2, 1, 1, 1;
1, 1, 1.2, 1, 1;
1, 1, 1.2, 0.8, 0.8];
loopNames = {'kDelta', 'kPhiDot', 'kPhi', 'kPhi'};
% the x and y indices for each plot
xy = [7, 15;
4, 12;
4, 12;
4, 12];
% where to get the x and y data from
xySource = {'outputs', 'outputsDot', 'outputs', 'outputs'};
% axis labels
xlabels = {'(a) $\delta$ (rad)',
'(b) $\dot{\phi}$ (rad/s)',
'(c) $\phi$ (rad)',
'(d) $\phi$ (rad)'};
ylabels = {'$\dot{\delta}$ (rad/s)',
'$\ddot{\phi}$ (rad/s$^2$)',
'$\dot{\phi}$ (rad/s)',
'$\dot{\phi}$ (rad/s)'};
% legend starts
legends = {'$k_\delta$ = ',
'$k_{\dot{\phi}}$ = ',
'$k_\phi$ = ',
'$k_\phi$ = '};
% how many decimals to show for each legend
floatSpec = {'%1.1f', '%1.3f', '%1.1f', '%1.1f'};
for i = 1:length(loopNames)
display(sprintf('Calculating phase portrait %d', i))
% adjust the gains and get the data
twentyPercent = generate_data('Benchmark', 5.0, ...
'gainMuls', gainChanges(i, :));
subplot(2, 2, i)
hold on
if i == 4
display('Phase portrait 4 comparison data.')
nominalData = generate_data('Benchmark', 5.0, ...
'gainMuls', [1, 1, 1, 0.8, 0.8]);
x = nominalData.(xySource{i})(:, xy(i, 1));
y = nominalData.(xySource{i})(:, xy(i, 2));
else
x = bikeData.(xySource{i})(:, xy(i, 1));
y = bikeData.(xySource{i})(:, xy(i, 2));
end
plot(x, y, 'k-', 'Linewidth', 1.0)
x = twentyPercent.(xySource{i})(:, xy(i, 1));
y = twentyPercent.(xySource{i})(:, xy(i, 2));
plot(x, y, 'k--', 'Linewidth', 1.0)
hold off
box on
axis equal
xlabel(xlabels{i}, 'Interpreter', 'Latex', 'Fontsize', 8)
ylabel(ylabels{i}, 'Interpreter', 'Latex', 'Fontsize', 8)
% make the legend
leg1 = sprintf(floatSpec{i}, bikeData.modelPar.(loopNames{i}));
leg2 = sprintf(floatSpec{i}, twentyPercent.modelPar.(loopNames{i}));
legend({[legends{i} leg1], [legends{i} leg2]} , ...
'Interpreter', 'Latex', ...
'Fontsize', 6)
end
plotAxes = findobj(gcf, 'Type', 'Axes');
set(plotAxes, 'Fontsize', 8)
% save the plot
filename = 'phasePortraits.eps';
print(gcf, ['plots' filesep filename], '-deps2c', '-loose')
fix_ps_linestyle(['plots' filesep filename])
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function eigenvalues(data, linestyles, colors)
global goldenRatio
figure()
figWidth = 5.0;
figHeight = figWidth / goldenRatio;
set(gcf, ...
'Color', [1, 1, 1], ...
'PaperOrientation', 'portrait', ...
'PaperUnits', 'inches', ...
'PaperPositionMode', 'manual', ...
'PaperPosition', [0, 0, figWidth, figHeight], ...
'PaperSize', [figWidth, figHeight])
bikes = fieldnames(data);
bikes(1) = [];
speeds = 0:0.1:10;
eVals = zeros(length(bikes), length(speeds), 11);
for i = 1:length(bikes)
% load the bicycle parameters
pathToParFile = ['parameters' filesep bikes{i} 'Par.txt'];
par = par_text_to_struct(pathToParFile);
str = 'Calculating eigenvalues for the %s bicycle and rider.';
display(sprintf(str, bikes{i}))
for j = 1:length(speeds)
% calculate the A, B, C, and D matrices of the bicycle model
[A, B, C, D] = whipple_pull_force_abcd(par, speeds(j));
eigenValues = eig(A);
eVals(i, j, :) = real(eig(A));
end
end
% reduce to the maximum values
zeroIndices = find(abs(eVals) <= 0.000001);
eVals(zeroIndices) = -100 * ones(size(zeroIndices));
maxEvals = max(eVals, [], 3);
lines = plot(speeds, maxEvals);
for i = 1:length(lines)
set(lines(i), ...
'Linestyle', linestyles{i}, ...
'Color', colors{i}, ...
'Linewidth', 1)
end
legend({'1', '2', '3', '4', '5', '6'})
xlabel('Speed (m/s)')
ylabel('Maximum real part of the eigenvalue (1/s)')
% set the tick marks differently
%set(gca, 'XTick', 0:0.5:10)
%set(gca, 'XTickLabel', {'0', '', '1', '', ...
%'2', '2.5', '3', '', ...
%'4', '', '5.0', '', ...
%'6', '', '7', '7.5', ...
%'8', '', '9', '', '10'})
%
% add some lines and labels for the speeds we looked at
hold on
maxLine = max(maxEvals, [], 1);
minLine = min(maxEvals, [], 1);
lines = zeros(4, 1);
speedInd = find(speeds == 2.5);
lines(1) = line([2.5, 2.5], ...
[minLine(speedInd) - 0.4, maxLine(speedInd) + 0.4]);
text(2, maxLine(speedInd) + 0.7, '2.5 m/s')
speedInd = find(speeds == 5.0);
lines(2) = line([5.0, 5.0], ...
[minLine(speedInd) - 0.4, maxLine(speedInd) + 0.4]);
text(4.5, maxLine(speedInd) + 0.7, '5.0 m/s')
speedInd = find(speeds == 7.5);
lines(3) = line([7.5, 7.5], ...
[minLine(speedInd) - 0.4, maxLine(speedInd) + 0.4]);
text(7, maxLine(speedInd) + 0.7, '7.5 m/s')
lines(4) = line([0, 10], [0, 0]);
hold off
set(lines, 'Color', 'k', 'Linewidth', 2)
% save the plot
filename = 'eigenvalues.eps';
print(gcf, ['plots' filesep filename], '-deps2c', '-loose')
fix_ps_linestyle(['plots' filesep filename])