using Convex using SCS using PyPlot PLOT_FIGURES = true # True to plot figures, false to suppress plots N = 96 # Number of periods in the day (so each interval 15 minutes) # Convenience variables for plotting fig_size = [14,3]; xtick_vals = [0, 8, 16, 24, 32, 40, 48, 56, 64, 72, 80, 88, 96]; xtick_labels = ["0:00", "2:00am", "4:00am", "5:00am", "8:00am", "10:00am", "12:00pm", "2:00pm", "4:00pm", "6:00pm", "8:00pm", "10:00pm", "12:00am"]; ############################################# # Price data generation - price values and intervals based off of PG&E Time Of Use plans ############################################# partial_peak_start = 34 # 08:30 peak_start = 48 # 12:00 peak_end = 72 # 18:00 (6:00pm) partial_peak_end = 86 # 21:30 (9:30pm) off_peak_inds = vcat(1:partial_peak_start, partial_peak_end+1:N); partial_peak_inds = vcat(partial_peak_start+1:peak_start, peak_end+1:partial_peak_end); peak_inds = vcat(peak_start+1:peak_end); # rates in $ / kWh off_peak_buy = 0.14; partial_peak_buy = 0.25; peak_buy = 0.45; # Rate cuts from buy prices to get sell prices off_peak_perc_cut = 0.20; partial_peak_perc_cut = 0.12; peak_perc_cut = 0.11; off_peak_sell = (1 - off_peak_perc_cut) * off_peak_buy; partial_peak_sell = (1 - partial_peak_perc_cut) * partial_peak_buy; peak_sell = (1 - peak_perc_cut) * peak_buy; # Combine the buy and sell prices into the price vectors R_buy = zeros(N); R_buy[off_peak_inds] .= off_peak_buy; R_buy[partial_peak_inds] .= partial_peak_buy; R_buy[peak_inds] .= peak_buy; R_sell = zeros(N); R_sell[off_peak_inds] .= off_peak_sell; R_sell[partial_peak_inds] .= partial_peak_sell; R_sell[peak_inds] .= peak_sell; # Plot the prices if PLOT_FIGURES figure(figsize=fig_size); plot(R_buy, label="Buy Price"); plot(R_sell, label="Sell Price"); legend(); ylabel("Price (\$ / kWh)"); title("Energy Prices (\$ / kWh)", fontsize=19); xticks(xtick_vals, xtick_labels ); end; ############################################# # Solar data generation ############################################# # Just something simple: a shifted cosine wave, squared to smooth edges, peak at noon # Some random +1 and +2 here and there to make it match python code exactly shift = (N/2); my_vec = ((1:N) .- shift .- 1) .* 2 .* pi ./ N p_pv = cos.(my_vec).^2; scale_factor = 35 p_pv = p_pv * scale_factor p_pv = max.(p_pv,0) p_pv[1:Int(shift/2)+1] .= 0 p_pv[N-Int(shift/2)+2:N] .= 0 # Plot it if PLOT_FIGURES figure(figsize=fig_size) plot(p_pv) title("PV Curve (kW)", fontsize=19) ylabel("Power (kW)") xticks(xtick_vals, xtick_labels ) end; ############################################# # Load Data Generation (using cvx) ############################################# # Fit a curve to some handpicked points and constrain the end points # to match and the derivative at the end to be the same at the beginning # points to fit to points = [ [1, 7], [11, 8], [21, 10], [29, 15], [37, 21], [46, 23], [53, 21], [57, 18], [61, 22.5], [67, 24.3], [71, 25], [74, 24], [84, 19], [96, 7], ] # Formulate an optimization problem that minimizes the error of # the fit while also minimizing the 2nd order difference of the function p_ld = Variable(N) obj_val = Variable(N + length(points)) prob = minimize(sum(obj_val)) # constr = [] prob.constraints += p_ld[1] == p_ld[N]; prob.constraints += (p_ld[2] - p_ld[1]) == (p_ld[N] - p_ld[N-1]); # Add a squared fitting error to the objective for i=1:length(points) prob.constraints += obj_val[i] >= square(p_ld[points[i][1]] - points[i][2]); end # Add the second order difference to the objective for i=2:N-1 prob.constraints += obj_val[i+length(points)] >= 100*square(p_ld[i+1] - 2*p_ld[i] + p_ld[i-1]); end prob.constraints += obj_val[1+length(points)] >= 100*square(p_ld[2] - 2*p_ld[1] + p_ld[N]); prob.constraints += obj_val[N+length(points)] >= 100*square(p_ld[1] - 2*p_ld[N] + p_ld[N-1]); solver = SCSSolver(verbose=0) solve!(prob, solver) p_ld = p_ld.value # Plot the curve if PLOT_FIGURES figure(figsize=fig_size) plot(p_ld) title("Load Curve (kW)", fontsize=19) ylabel("Power (kW)") xticks(xtick_vals, xtick_labels ) end; # # For reference, plot the net load curve: load - solar # # Aside: the shape of this curve is referred to as "the duck curve" # # by some people in energy (since the shape looks like a duck on profile), and # # is a result of renewable generation displacing load demand in the day. # if PLOT_FIGURES # figure(figsize=fig_size) # plot(p_ld - p_pv) # axhline(0, color="black", linestyle="--") # title("Load minus Solar (kW)", fontsize=19) # ylabel("Power (kW)") # xticks(xtick_vals, xtick_labels ) # show() # end ############################################# # Battery and Grid Line Constraint Values ############################################# # Max charge and discharge rates D = 10; # Max discharge rate (kW) C = 8; # Max charge rate (kW) Q = 27; # Max energy (kWh) # Final list of values generated: # N (scalar): number of intervals we split the day into # R_buy (vector, $/kWh): prices one can buy energy at from grid in given interval # R_sell (vector, $/kWh): prices one can sell energy at to grid in given interval # p_pv (vector, kW): power generated by solar # p_ld (vector, kW): power demands of load # D (scalar, kW): max discharge rate of battery # C (scalar, kW): max charge rate of battery # Q (scalar, kWh): max energy of battery