file cleanup and reorg

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

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])

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);