From 94ba01895bcc8e07c40691896849daf8976079cc Mon Sep 17 00:00:00 2001 From: Kevin Moy <68761709+kevinrussellmoy@users.noreply.github.com> Date: Thu, 8 Dec 2022 17:49:50 -0800 Subject: [PATCH] file cleanup and reorg --- genset_data.m => data/genset_data.m | 0 load_data.m => data/load_data.m | 0 pv_data.m => data/pv_data.m | 0 pmp_control_ECM.m | 258 +----------------- .../fuelcons_degrad_output.mat | Bin .../soc_fuelcons_degrad_20210920.mat | Bin .../soc_fuelcons_degrad_20211007.mat | Bin .../soc_fuelcons_degrad_toy_20210920.mat | Bin soc_fuelcons_degrad.m | 82 ------ soc_fuelcons_degrad_toy.m | 9 - 10 files changed, 1 insertion(+), 348 deletions(-) rename genset_data.m => data/genset_data.m (100%) rename load_data.m => data/load_data.m (100%) rename pv_data.m => data/pv_data.m (100%) rename fuelcons_degrad_output.mat => results/fuelcons_degrad_output.mat (100%) rename soc_fuelcons_degrad_20210920.mat => results/soc_fuelcons_degrad_20210920.mat (100%) rename soc_fuelcons_degrad_20211007.mat => results/soc_fuelcons_degrad_20211007.mat (100%) rename soc_fuelcons_degrad_toy_20210920.mat => results/soc_fuelcons_degrad_toy_20210920.mat (100%) diff --git a/genset_data.m b/data/genset_data.m similarity index 100% rename from genset_data.m rename to data/genset_data.m diff --git a/load_data.m b/data/load_data.m similarity index 100% rename from load_data.m rename to data/load_data.m diff --git a/pv_data.m b/data/pv_data.m similarity index 100% rename from pv_data.m rename to data/pv_data.m diff --git a/pmp_control_ECM.m b/pmp_control_ECM.m index 1faf4d8..8132b40 100644 --- a/pmp_control_ECM.m +++ b/pmp_control_ECM.m @@ -164,7 +164,7 @@ % control given costate % Costate -lambda_inits = linspace(-150, -50, 41); +lambda_inits = linspace(-150, -50, 21); for i = 1:length(lambda_inits) [u_opt, x, lambda] = mgpmpecm(length(ld1), pv1, ld1, x_max, x_min, LIB_INV_SIZE_KW, LIB_EFF_CHG, ... @@ -269,259 +269,3 @@ legend(legendStrings, 'location', 'eastoutside') set(gca, "FontSize", 20) - -%% SCRATCH - -% yr_dt_sec = (YEAR_START:seconds(1):datetime(2021,12,31,23,59,59))'; - -% week_len = 4*24*7; -% day_len = 4*24; -% % week_start = 0 -% week_start = 30 * week_len; % try 25? -% week_end = week_start + week_len-1; -% -% day_start = week_start; -% day_end = day_start + day_len-1; - -% % One week: -% ld_wk1 = ld(week_start:week_end); -% pv_wk1 = pv(week_start:week_end); -% wk1_dt = yr_dt(week_start:week_end); -% -% % One day (START WITH THIS): -% ld1 = ld(day_start:day_end); -% pv1 = pv(day_start:day_end); -% d1_dt = yr_dt(day_start:day_end); - -% % Convert to seconds -% ld_sec = repelem(ld1, 3600*h); -% pv_sec = repelem(pv1, 3600*h); -% d1_dt_sec = yr_dt_sec(30*24*7*3600+1:30*24*7*3600+(3600*24)); - -% %% 1st simple case: Fixed energy, power of ESS; fixed week, fixed costate; fixed starting energy -% % Initial value of costate -% lambda_init = 10; -% -% [u_opt, x, lambda] = mgpmp(length(ld_sec), pv_sec, ld_sec, x_max, x_min, LIB_INV_SIZE_KW, LIB_EFF_CHG, ... -% socdot, P_batt_range_socdot, SOC_range_socdot, ... -% dfdx, P_batt_range_dfdx, SOC_range_dfdx, ... -% h, lambda_init); -% -% % Get resulting curtailed PV, DG generation, final energy at optimal -% % control given costate -% dg1 = max(ld_sec-pv_sec-u_opt,0); -% pv_curt = max(pv_sec-ld_sec-u_opt,0); -% x_f = x(end); -% -% % Plot Optimal Dispatch!! -% hFig = figure(3); -% set(hFig, 'Position', [100 100 2000 2200]) -% subplot(2,2,1) -% hold on -% plot(d1_dt_sec, ld_sec, 'LineWidth', 2) -% plot(d1_dt_sec, pv_sec, 'LineWidth', 2) -% plot(d1_dt_sec, u_opt, 'LineWidth', 2) -% plot(d1_dt_sec, dg1, 'LineWidth', 2) -% ylabel("Power [kWAC]") -% legend("Load", "PV", "Optimal ESS Dispatch", "Diesel Genset", "location", "best") -% set(gca, "FontSize", 20) -% subplot(2,2,2) -% plot(yr_dt_sec(30*24*7*3600+1:30*24*7*3600+(3600*24)+1), x, 'LineWidth', 2) -% ylim([0 1]) -% ylabel("SOC [-]") -% set(gca, "FontSize", 20) -% subplot(2,2,3) -% plot(yr_dt_sec(30*24*7*3600+1:30*24*7*3600+(3600*24)+1), lambda, 'LineWidth', 2) -% ylabel("\lambda [L]") -% set(gca, "FontSize", 20) -% subplot(2,2,4) -% plot(d1_dt_sec, cumsum(genset_model(dg1))*h/3600, 'LineWidth', 2) -% ylabel('Cumulative Fuel Consumption [L]') -% set(gca, "FontSize", 20) - -%% Function to calculate optimal control, final state, and state vector -% % TODO: DELETE THIS HERE, USE mgpmpecm.m INSTEAD -% function [u_opt, x, lambda] = mgpmp(outage_len, pv_outage, ld_outage, ... -% x_max, x_min, inv_rating, inv_eff, ... -% socdot, P_batt_range_socdot, SOC_range_socdot, ... -% dfdx, P_batt_range_dfdx, SOC_range_dfdx, ... -% h, lambda_init) -% % Function to compute optimal control strategy based on Pontryagin -% % Minimization Principle (PMP) -% % Inputs: -% % outage_len = length of outage period (e.g. one day) in 15 min. intervals -% % pv_outage = pv profile during outage period, as kW -% % ld_outage = load profile during outage period, as kW -% % x_max = upper bound of state -% % x_min = lower bound of state -% % inv_rating = kWAC rating of LIB inverter -% % inv_eff = AC/DC conversion efficiency of LIB inverter -% % socdot = map from (pack power, SOC) --> change in SOC -% % P_batt_range_socdot = corresponding pack power range for socdot map -% % SOC_range_socdot = corresponding SOC range for socdot map -% % dfdx = map from (pack power, SOC) --> dSOC_dot/dSOC -% % P_batt_range_dfdx = corresponding pack power range for dfdx map -% % SOC_range_dfdx = corresponding SOC range for dfdx map -% % lambda = initial costate -% % h = fraction of the hour -% % Outputs: -% % u_opt = optimal control -% % x = state as a function of time from t=1:day_len+1 -% % lambda = costate -% x = zeros(outage_len+1, 1); -% u_opt = zeros(outage_len, 1); -% lambda = zeros(outage_len+1, 1); -% -% % TODO: Change this assumption -% x(1) = x_max; -% -% lambda(1) = lambda_init; -% -% % Set penalty (K) for violating constraints -% K = 200; % TODO: Sensitivity analysis on this -% -% % for t = 1:3600*5 -% for t = 1:outage_len -% % disp('Time') -% % disp(t) -% % First, work in all AC power -% % Determine control space at time t based on PV, load power -% % assuming that ESS can only charge off of PV, discharge to load -% % But bounded by inverter rating -% u_min_t = -min(pv_outage(t), inv_rating); -% u_max_t = min(ld_outage(t), inv_rating); -% % u_min_t = -pv_outage(t); -% % u_max_t = ld_outage(t); -% u_range_t = sort([linspace(u_min_t,u_max_t) 0]); % include 0 power as an option -% % disp(u_range_t) -% -% % Calculate optimal control with control space at time t -% H_t = zeros(length(u_range_t), 1); -% wE_t = zeros(length(u_range_t), 1); -% -% for v = 1:length(u_range_t) -% % Put rules (~ environment) for each dispatch here -% % TODO: Break this out into its own function (or file/method??) -% % Initialize PV, load, and ESS energy resources; DG power, penalty -% pv_t = pv_outage(t); -% l_t = ld_outage(t); -% x_t = x(t); -% dg_t = 0; -% -% u_t = u_range_t(v); -% -% % First, charge PV off of ESS if u is negative -% pv_t = pv_t - (-min(0,u_t)); -% -% % Dispatch portion of load from ESS if u is positive -% l_t = l_t - (max(0,u_t)); -% -% % Balance the load and PV, with excess PV curtailed or excess load -% % supplied by DG -% if pv_t > l_t -% % Subtract off load from PV and curtail the rest -% % DG is not used here -% elseif pv_t < l_t -% % Subtract off PV from load and supply remainder load with DG -% l_t = l_t - pv_t; -% dg_t = l_t; -% else -% % Load and PV are in perfect balance! How?? -% % DG is not used here -% end -% % disp(dg_t) -% % Calculate penalty contribution -% % First, convert AC power to DC power -% if u_t >= 0 -% u_DC_t = u_t/inv_eff; -% else -% u_DC_t = u_t * inv_eff; -% end -% -% % Use socdot map to find SOC at time step t+1 -% socdot_t = interp2(P_batt_range_socdot, SOC_range_socdot, socdot, u_DC_t, x_t, 'spline'); -% -% x_t1 = x_t + socdot_t*h; -% % disp(x_t1) -% % x_t1 = x_t + socdot_t; -% if x_t1 > x_max -% w = K; -% elseif x_t1 < x_min -% w = -K; -% elseif x_min <= x_t1 && x_t1 <= x_max -% w = 0; -% % else -% % w = K; -% end -% % disp(w) -% wE_t(v) = w; -% % disp(w) -% % Calculate fuel consumpt ion from function -% fuel_cons_L = genset_model(dg_t); -% % fuel_cons_L = genset_model(dg_t) * h / 3600; -% -% % Put resulting Hamiltonian here -% % H_t(v) = fuel_cons_L + (lambda(t) + w) * u_t; -% H_t(v) = fuel_cons_L + (lambda(t) + w) * socdot_t; -% % disp(H_t(v)) -% end -% % % Plot -% % figure(2) -% % subplot(2,1,1) -% % plot(u_range_t,H_t) -% % ylabel("Hamiltonian value") -% % xlabel("Power, kW") -% % subplot(2,1,2) -% % plot(u_range_t, wE_t, 'o') -% % ylabel("Penalty value") -% % xlabel("Power, kW") -% -% % Get index of H_t where it is minimized; -% [~, min_ind] = min(H_t); -% % disp(min(H_t)) -% % disp(min_ind) -% % Find the corresponding LIB AC power value as the optimal control -% u_opt_t = u_range_t(min_ind); -% -% % Store optimal LIB AC power control -% u_opt(t) = u_opt_t; -% -% % First, convert AC power to DC power -% if u_opt_t >= 0 -% u_DC_opt_t = u_opt_t/inv_eff; -% else -% u_DC_opt_t = u_opt_t * inv_eff; -% end -% -% % Use socdot map to find SOC at time step t+1 -% socdot_t = interp2(P_batt_range_socdot, SOC_range_socdot, socdot, u_DC_opt_t, x(t), 'spline'); -% -% % Update the SOC for the next time step -% % x(t+1) = x(t) + socdot_t; -% x(t+1) = x(t) + socdot_t*h; -% % disp(x(t+1)) -% % Update the costate -% % 1) Calculate penalty -% if x(t) > x_max -% w = K; -% elseif x(t) < x_min -% w = -K; -% elseif x_min <= x(t) && x(t) <= x_max -% w = 0; -% % else -% % w = K; -% end -% -% % 2) Compute dfdx -% dfdx_t = -interp2(P_batt_range_dfdx, SOC_range_dfdx, dfdx, u_DC_opt_t, x(t), 'spline'); -% -% % 3) Compute lambda_dot -% lambda_dot = -(lambda(t) + w)*dfdx_t; -% -% % 4) Evolve costate -% lambda(t+1) = lambda(t) + lambda_dot*h; -% end -% -% -% % function end -% end diff --git a/fuelcons_degrad_output.mat b/results/fuelcons_degrad_output.mat similarity index 100% rename from fuelcons_degrad_output.mat rename to results/fuelcons_degrad_output.mat diff --git a/soc_fuelcons_degrad_20210920.mat b/results/soc_fuelcons_degrad_20210920.mat similarity index 100% rename from soc_fuelcons_degrad_20210920.mat rename to results/soc_fuelcons_degrad_20210920.mat diff --git a/soc_fuelcons_degrad_20211007.mat b/results/soc_fuelcons_degrad_20211007.mat similarity index 100% rename from soc_fuelcons_degrad_20211007.mat rename to results/soc_fuelcons_degrad_20211007.mat diff --git a/soc_fuelcons_degrad_toy_20210920.mat b/results/soc_fuelcons_degrad_toy_20210920.mat similarity index 100% rename from soc_fuelcons_degrad_toy_20210920.mat rename to results/soc_fuelcons_degrad_toy_20210920.mat diff --git a/soc_fuelcons_degrad.m b/soc_fuelcons_degrad.m index 491ff23..cb32561 100644 --- a/soc_fuelcons_degrad.m +++ b/soc_fuelcons_degrad.m @@ -33,15 +33,6 @@ outage_lens(i) = length(ld(t_ind(1):t_ind(2))); end -% % Select load and PV data -% -% ld1 = ld(t_ind(1):t_ind(2)); -% pv1 = pv(t_ind(1):t_ind(2)); -% dt1 = yr_dt(t_ind(1):t_ind(2)); -% -% % Create variable for outage length (v useful) -% outage_len = length(ld1); - %% SOCs to try SOCs = linspace(x_min, x_max, 27); @@ -115,79 +106,6 @@ cap_pcts{k} = cap_losses/Q_nom; end -% %% Search for optimal initial costate values -% lambda_min = -150; -% lambda_max = -75; -% disp('Finding optimal initial costate values for each value of SOC:') -% for i = 1:length(SOCs) -% disp(['Iteration ', num2str(i)]) -% [s_init_vec, SOC_end_vec] = init_costate_search(length(ld1), pv1, ld1, x_max, x_min, ... -% LIB_INV_SIZE_KW, LIB_EFF_CHG, ... -% socdot, P_batt_range_socdot, SOC_range_socdot, ... -% dfdx, P_batt_range_dfdx, SOC_range_dfdx, ... -% h, lambda_min, lambda_max, SOCs(i)); -% -% SOC_fins(i) = SOC_end_vec(end); -% lambda_inits(i) = s_init_vec(end); -% disp(['SOC ', num2str(SOCs(i)), ' Optimal costate value ', num2str(lambda_inits(i)), ' found']) -% end -% disp('Done finding optimal initial costate values!') -% -% %% Plot -% hFig = figure(1); -% % set(hFig, 'Position', [100 100 920 500]) -% set(hFig, 'Position', [100 100 600 500]) -% hold on -% plot(SOC_fins, lambda_inits, 'marker','^', 'MarkerSize',15) -% xlabel("Initial Costate Value [L/hr]") -% ylabel("Final SOC [-]") -% set(gca, "FontSize", 28) -% -% %% Find fuel consumption \dot{m_f}(t), LIB power output u(t), and LIB SOC(t) at each optimal initial costate value -% % TODO: Combine with above so we are not doing redundant calculations -% disp('Finding optimal fuel consumption, SOC, LIB power from optimal intial costate values:') -% for i = 1:length(lambda_inits) -% disp(['Iteration ', num2str(i)]) -% [u_opt, x, lambda] = mgpmpecm(length(ld1), pv1, ld1, x_max, x_min, LIB_INV_SIZE_KW, LIB_EFF_CHG, ... -% socdot, P_batt_range_socdot, SOC_range_socdot, ... -% dfdx, P_batt_range_dfdx, SOC_range_dfdx, ... -% h, lambda_inits(i)); -% dg1 = max(ld1-pv1-u_opt,0); -% pv_curt = max(pv1-ld1-u_opt,0); -% x_f = x(end); -% -% u_opts(:,i) = u_opt; -% SOC_opts(:,i) = x; -% fuel_consumps(i) = sum(genset_model(dg1))*h; -% -% end -% disp('Done with that.') -% -% % Compute baseline fuel consumption -% L_base = sum(genset_model(max(ld1-pv1,0)))*h; -% -% %% -% disp('Finding degradation from optimal SOC trajectories:') -% for i = 1:length(SOCs) -% disp(['Iteration ', num2str(i)]) -% -% %TODO WHY IS THIS IMAGINARY?? -% [cap_loss] = emp_deg_model(SOC_opts(:,i), Q_nom, h); -% cap_losses(:,i) = cap_loss; -% end -% disp('Done with that.') -% %% Plot -% -% hFig = figure(2); -% % set(hFig, 'Position', [100 100 920 500]) -% set(hFig, 'Position', [100 100 600 500]) -% hold on -% plot(SOCs, cap_losses, 'marker','.', 'MarkerSize',15, 'LineWidth', 2) -% xlabel("Final SOC [-]") -% xlim([0.2 0.85]) -% ylabel("Capacity Loss [Ah]") -% set(gca, "FontSize", 28) - %% Plot Pareto Front %indices to plot given restrictions on Andrea's model %(SOC_min = [0.2,0.45]) diff --git a/soc_fuelcons_degrad_toy.m b/soc_fuelcons_degrad_toy.m index 31447f3..ad3770d 100644 --- a/soc_fuelcons_degrad_toy.m +++ b/soc_fuelcons_degrad_toy.m @@ -33,15 +33,6 @@ outage_lens(i) = length(ld(t_ind(1):t_ind(2))); end -% % Select load and PV data -% -% ld1 = ld(t_ind(1):t_ind(2)); -% pv1 = pv(t_ind(1):t_ind(2)); -% dt1 = yr_dt(t_ind(1):t_ind(2)); -% -% % Create variable for outage length (v useful) -% outage_len = length(ld1); - %% SOCs to try SOCs = linspace(x_min, x_max, 27);