Merge pull request #8 from itchono/asc2024

ASC2024

Author: itchono

.gitignore

@@ -13,4 +13,7 @@ resources/**
 *.txt
 
 # Animations
-anim
+anim
+
+# Workspaces
+*.mat

postprocessing/plotting/imseq_vecex.m

@@ -0,0 +1,81 @@
+function imseq_vecex(y, cfg, nfmax, view_dim)
+% plot for vector export, based on maximum # of vector lines and frames
+% Expect y to be (6, N)
+% no need to interpolate vectors unless you have VERY big skips in L
+% Reference: https://www.mathworks.com/help/matlab/ref/getframe.html
+% Splits by orbit number
+
+[~, ~, ~, ~, ~, L] = unpack_mee(y);
+
+if nargin < 3
+    nfmax = 20;
+    view_dim = 3;
+end
+
+%% Data processing
+ind_orbits = find(diff(mod(L, 2*pi)) < 0);
+num_orbits = length(ind_orbits);
+
+%% Initial plot stuff
+L_sample = linspace(0, 2*pi, 60);
+plot_sphere(0, 0, 0, 6378e3, [0.3010, 0.7450, 0.9330]);
+hold on
+
+cm = colormap("turbo");
+th = title(sprintf("Orbit 1 of %d", num_orbits));
+
+% plot initial trace
+cart_sample = mee2cartesian([repmat(y(1:5, 1), 1, numel(L_sample)); L_sample]);
+plot3(cart_sample(1, :), cart_sample(2, :), cart_sample(3, :), "Color", cm(1, :), "DisplayName", "Initial Orbit", "LineWidth", 1, "LineStyle", "--");
+
+% plot main orbit
+main_orbit = plot3(cart_sample(1, :), cart_sample(2, :), cart_sample(3, :), "Color", "black", "DisplayName", "Current Orbit", "LineWidth", 1);
+
+% plot final trace
+cart_sample = mee2cartesian([repmat(y(1:5, end), 1, numel(L_sample)); L_sample]);
+plot3(cart_sample(1, :), cart_sample(2, :), cart_sample(3, :), "Color", cm(end, :), "DisplayName", "Final Orbit", "LineWidth", 1, "LineStyle", "--");
+
+% colorbar
+colorbar;
+clim([1,num_orbits]); 
+
+axis equal
+view(view_dim)
+
+%% Initialize path
+mkdir("anim")
+mkdir("anim", cfg.casename)
+
+%% Animation
+% spacing; never plot more than 20 frames
+stride = max(1, ceil(num_orbits/nfmax));
+
+% ensure we always plot the last orbit
+ii = 1; % counter
+for j = [1:stride:num_orbits, num_orbits]
+    % Update plots
+    idx = ind_orbits(j);
+    cart_sample = mee2cartesian([repmat(y(1:5, idx), 1, numel(L_sample)); L_sample]);
+    main_orbit.XData = cart_sample(1, :);
+    main_orbit.YData = cart_sample(2, :);
+    main_orbit.ZData = cart_sample(3, :);
+    
+    % Add shadow of previous orbits
+    colour_idx = ceil(j/num_orbits*length(cm));
+    plot3(cart_sample(1, :), cart_sample(2, :), cart_sample(3, :), "Color", cm(colour_idx, :), "LineWidth", 1);
+    
+    th.String = sprintf("Orbit %d of %d", j, num_orbits);
+    
+    % Draw and get frame
+    drawnow
+
+    % https://www.mathworks.com/matlabcentral/answers/363832-some-figures-not-saving-as-vector-graphics-svg
+    % https://www.mathworks.com/matlabcentral/answers/471164-print-pdf-to-a-specific-size
+    hf = gcf;
+    hf.Units = 'centimeters';        % set figure units to cm
+    hf.PaperUnits = 'centimeters';   % set pdf printing paper units to cm
+    hf.PaperSize = hf.Position(3:4);  % assign to the pdf printing paper the size of the figure
+    print(hf,'-vector','-dpdf',fullfile("anim", cfg.casename, sprintf("anim-%d.pdf", ii))) % pdf
+    ii = ii + 1; % JANK but works
+end
+end

postprocessing/plotting/plot_orbit_vecex.m

@@ -0,0 +1,53 @@
+% plot for vector export, maximum # of lines etc.
+
+function plot_orbit_vecex(y)
+% Expect y to be (6, N)
+% Splits by orbit number so that maximum number of lines is under 300
+
+[~, ~, ~, ~, ~, L] = unpack_mee(y);
+
+%% Data processing
+ind_orbits = find(diff(mod(L, 2*pi)) < 0);
+num_orbits = length(ind_orbits);
+
+%% Initial plot stuff
+L_sample = linspace(0, 2*pi, 60);
+plot_sphere(0, 0, 0, 6378e3, [0.3010, 0.7450, 0.9330]);
+hold on
+
+cm = colormap("turbo");
+
+NMAX = 200;
+
+% spacing
+stride = max(1, ceil(num_orbits/NMAX));
+
+% ensure we always plot the last orbit
+for j = [1:stride:num_orbits, num_orbits]
+    % Update plots
+    idx = ind_orbits(j);
+    cart_sample = mee2cartesian([repmat(y(1:5, idx), 1, numel(L_sample)); L_sample]);
+    
+    % Add shadow of previous orbits
+    colour_idx = ceil(j/num_orbits*length(cm));
+    plot3(cart_sample(1, :), cart_sample(2, :), cart_sample(3, :), "Color", cm(colour_idx, :), "LineWidth", 0.5);
+end
+
+title(sprintf("Earth Inertial Coordinates (Showing Every %d Orbits)", stride));
+xlabel("x (m)")
+ylabel("y (m)")
+zlabel("z (m)")
+cbar = colorbar;
+
+% wrangle the discrete colorbar tick alignment
+cbscalein = cbar.Limits;
+cbscaleout = [0 num_orbits];
+nt = 8;
+ticks = linspace(cbscaleout(1),cbscaleout(2), nt);
+cbar.Ticks = diff(cbscalein)*(ticks-cbscaleout(1))/diff(cbscaleout) + cbscalein(1);
+cbar.TickLabels = round(linspace(0, num_orbits, nt)');
+
+ylabel(cbar, "Orbit Number")
+axis equal;
+view(3)
+end

postprocessing/postprocess.m

@@ -6,20 +6,26 @@
 plot_elements_mee(y, t, cfg.y_target);
 exportgraphics(hf1, fullfile("outputs", casename, 'orbital_elements.pdf'), 'ContentType', 'vector')
 
-hf2 = figure;
-plot_elements_ke(y, t, cfg.y_target);
-exportgraphics(hf2, fullfile("outputs", casename, 'orbital_elements_ke.pdf'), 'ContentType', 'vector')
+% KE plot is not very critical
+% hf2 = figure;
+% plot_elements_ke(y, t, cfg.y_target);
+% exportgraphics(hf2, fullfile("outputs", casename, 'orbital_elements_ke.pdf'), 'ContentType', 'vector')
 
-% Split up plots
-hf3 = figure;
-plot_steering_history(y, t, cfg);
-exportgraphics(hf3, fullfile("outputs", casename, 'steering_history.pdf'), 'ContentType', 'vector')
+% old plot - nothing very useful came of it for production purposes
+% hf3 = figure;
+% plot_steering_history(y, t, cfg);
+% exportgraphics(hf3, fullfile("outputs", casename, 'steering_history.pdf'), 'ContentType', 'vector')
 
-hf4 = figure;
-[y_interp, ~] = interp_mee(y, t, 100);
-plot_orbit_mee(y_interp);
-exportgraphics(hf4, fullfile("outputs", casename, 'trajectory_plot.png'), ...
-    'Resolution', 300)
+% bitmap plot superceded by vector output
+% hf4 = figure;
+% [y_interp, ~] = interp_mee(y, t, 100);
+% plot_orbit_mee(y_interp);
+% exportgraphics(hf4, fullfile("outputs", casename, 'trajectory_plot.png'), ...
+%     'Resolution', 300)
+
+hf5 = figure;
+plot_orbit_vecex(y);
+exportgraphics(hf5, fullfile("outputs", casename, 'trajectory_plot.pdf'), 'ContentType', 'vector')
 
 %% Video
 anim_osculating_mee(y, fullfile("outputs", casename, 'trajectory_animation.mp4'));

scripts/benchmark_transfer.m

@@ -3,7 +3,7 @@
 % Supposed to be a worst-case for the guidance law
 % As of the latest version, we get a performance of 874 revolutions, and
 % about 6.7 km/s of delta-v expenditure using the unpenalized law
-% to within a tolerance of 1e-3; somewhat sesitive to NDF angles
+% to within a tolerance of 1e-3; somewhat sensitive to NDF angles
 
 %% Problem Definition
 % Create a struct for neatness
@@ -14,7 +14,7 @@
 cfg.solver = @ode89;
 cfg.t_span = [0, 1e8];
 cfg.options = odeset('RelTol', 1e-4, "Stats", "on", "MaxStep", 1e4);
-cfg.tol = 1e-3;
+cfg.tol = 5e-3;
 cfg.guidance_weights = [1; 1; 1; 1; 1];
 cfg.penalty_param = 1;
 cfg.min_pe = 10000e3;

scripts/benchmark_transfer_j2.m

@@ -1,5 +1,5 @@
 %% Description
-% J2 version of benchmark transfer
+% Determine effects of J2 on benchmark transfer
 
 %% Problem Definition
 % Create a struct for neatness
@@ -10,7 +10,7 @@
 cfg.solver = @ode89;
 cfg.t_span = [0, 1e8];
 cfg.options = odeset('RelTol', 1e-4, "Stats", "on", "MaxStep", 1e4);
-cfg.tol = 1e-3;
+cfg.tol = 5e-3;
 cfg.guidance_weights = [1; 1; 1; 1; 1];
 cfg.penalty_param = 1;
 cfg.min_pe = 10000e3;

scripts/benchmark_transfer_optim.m renamed to scripts/benchmark_transfer_optim_100.m

@@ -11,11 +11,11 @@
 cfg.t_span = [0, 1e8];
 cfg.options = odeset('RelTol', 1e-4, "Stats", "on", "MaxStep", 1e4);
 cfg.tol = 5e-3;
-cfg.guidance_weights = [1.774e+00; 5.149e-01; 3.327e-01; 9.925e+00; 5.317e-01];
+cfg.guidance_weights = [ 6.810e-01;  3.343e-01;  4.072e-01;  9.965e+00; 4.360e-01];
 cfg.penalty_param = 1;
 cfg.min_pe = 10000e3;
 cfg.penalty_weight = 0;
-cfg.kappa = deg2rad(64);
+cfg.kappa = deg2rad(61.3);
 cfg.dynamics = "mee";
 cfg.j2 = false;
 

scripts/benchmark_transfer_optim_20.m

@@ -0,0 +1,24 @@
+%% Description
+% Version of benchmark transfer with optimized weights
+
+%% Problem Definition
+% Create a struct for neatness
+cfg.y0 = [20000e3; 0.5; -0.2; 0.5; 0; 0];
+cfg.y_target = [25000e3; 0.2; 0.5; 0; 0.3];
+cfg.propulsion_model = @sail_thrust;
+cfg.steering_law = @quail;
+cfg.solver = @ode89;
+cfg.t_span = [0, 1e9];
+cfg.options = odeset('RelTol', 1e-4, "Stats", "on", "MaxStep", 1e4);
+cfg.tol = 5e-3;
+cfg.guidance_weights = [2.06270441;  0.95005485;  4.44929224;  6.58042015;  7.3647472 ];
+cfg.penalty_param = 1;
+cfg.min_pe = 10000e3;
+cfg.penalty_weight = 0;
+cfg.kappa = deg2rad(56.32);
+cfg.dynamics = "mee";
+cfg.j2 = false;
+
+%% Run
+[~, cfg.casename, ~] = fileparts(mfilename);
+[y, t, dv] = run_mission(cfg);

scripts/benchmark_transfer_optim_j2.m renamed to scripts/benchmark_transfer_optim_200.m

@@ -11,13 +11,14 @@
 cfg.t_span = [0, 1e8];
 cfg.options = odeset('RelTol', 1e-4, "Stats", "on", "MaxStep", 1e4);
 cfg.tol = 5e-3;
-cfg.guidance_weights = [1.774e+00; 5.149e-01; 3.327e-01; 9.925e+00; 5.317e-01];
+cfg.guidance_weights = [ 1.485e+00;  3.472e+00;  1.973e-01;  7.157e+00;
+        	9.814e+00];
 cfg.penalty_param = 1;
 cfg.min_pe = 10000e3;
 cfg.penalty_weight = 0;
-cfg.kappa = deg2rad(64);
+cfg.kappa = deg2rad(59.83);
 cfg.dynamics = "mee";
-cfg.j2 = true;
+cfg.j2 = false;
 
 %% Run
 [~, cfg.casename, ~] = fileparts(mfilename);

scripts/benchmark_transfer_optim_500.m

@@ -0,0 +1,25 @@
+%% Description
+% Version of benchmark transfer with optimized weights
+
+%% Problem Definition
+% Create a struct for neatness
+cfg.y0 = [20000e3; 0.5; -0.2; 0.5; 0; 0];
+cfg.y_target = [25000e3; 0.2; 0.5; 0; 0.3];
+cfg.propulsion_model = @sail_thrust;
+cfg.steering_law = @quail;
+cfg.solver = @ode89;
+cfg.t_span = [0, 1e8];
+cfg.options = odeset('RelTol', 1e-4, "Stats", "on", "MaxStep", 1e4);
+cfg.tol = 1e-2;
+cfg.guidance_weights = [ 6.827e-01;  9.304e-01;  1.111e+00;  3.386e+00;
+        	8.306e+00 ];
+cfg.penalty_param = 1;
+cfg.min_pe = 10000e3;
+cfg.penalty_weight = 0;
+cfg.kappa = deg2rad(62.51);
+cfg.dynamics = "mee";
+cfg.j2 = false;
+
+%% Run
+[~, cfg.casename, ~] = fileparts(mfilename);
+[y, t, dv] = run_mission(cfg);