-
Notifications
You must be signed in to change notification settings - Fork 2
/
generate_all_figures.m
153 lines (140 loc) · 6.51 KB
/
generate_all_figures.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
%--------------------------------------------------------------%
% File: generate_all_figures.m (script)
% Author: Miguel Dias
% Date 25/08/16
% v1.0
% Description: Analysis of the solution found by the shortest_path solver:
% Plots the figures of tariff, demand and solution for all 36 cases.
% Exports the figures to eps format so they can be compiled into a single
% document using LaTex.
%--------------------------------------------------------------%
%Startup:
% The variable decided_costs takes too much room to be saved normally. If
% it already exists in the workspace, don't delete it. If it does not
% exist, load it and convert it to double.
close all;
clc;
if (exist('decided_costs', 'var'))
clearvars('-except', 'decided_costs');
load('..\Data\graph_data_all_days.mat', '-regexp', '[^decided_costs]');
else
clearvars;
load('..\Data\graph_data_all_days.mat');
decided_costs = double(decided_costs);
end
% d_index; 1-3-> Large; 4-6->FSR; 7-9->Small restaurant; % 10-12->Residential
%d_index=mod(index,12)+[mod(index,12)==0]*12;
end_time = 24; %end time in h
dt = 15; %in s
n_lines = 3600 / dt;
T = end_time * n_lines; %number of time steps
t_plot = (1:n_lines * end_time)' ./ n_lines;
savepath = '..\Data\Case_Study_Images\';
save_fig = 1;
%%
joule2kWh = 1 / 3.6e6; %1kWh=3.6e6J
fuel_index = [1 * ones(12, 1); 2 * ones(12, 1); 3 * ones(12, 1)];
%% Initialize all vectors for the loop
electric_energy_kWh = zeros(1, numel(fuel_index));
heat_energy_kWh = zeros(1, numel(fuel_index));
base_electricity_charge = zeros(1, numel(fuel_index));
base_heat_charge = zeros(1, numel(fuel_index));
total_charge = zeros(1, numel(fuel_index));
power_MGT = zeros(T, numel(fuel_index));
heat_MGT = zeros(T, numel(fuel_index));
mdot_MGT = zeros(T, numel(fuel_index));
new_demand = zeros(T, numel(fuel_index));
path_cost = zeros(1, numel(fuel_index));
bought_elec = zeros(1, numel(fuel_index));
sold_energy = zeros(1, numel(fuel_index));
bought_fuel = zeros(1, numel(fuel_index));
bought_heat = zeros(1, numel(fuel_index));
MGT_cost = zeros(1, numel(fuel_index));
savings = zeros(1, numel(fuel_index));
FC = zeros(1, numel(fuel_index));
MGT_PDC = zeros(1, numel(fuel_index));
ut_PDC = zeros(1, numel(fuel_index));
MGT_IDC = zeros(1, numel(fuel_index));
ut_IDC = zeros(1, numel(fuel_index));
%%
FS = 20;
for i = 1:numel(fuel_index)
d_index = mod(i, 12) + (mod(i, 12) == 0) * 12;
%Compute base costs
electric_energy_kWh(:, i) = sum(power_demand(:, d_index).*dt) * joule2kWh; %sum(w*s)= energy in kWh
heat_energy_kWh(:, i) = sum(heat_demand(:, d_index)*dt) * joule2kWh;
base_electricity_charge(:, i) = sum((power_demand(:, d_index) .* dt).*joule2kWh.*elec_tariff(:, d_index));
base_heat_charge(:, i) = sum((heat_demand(:, d_index) .* dt).*joule2kWh.*heat_tariff(fuel_index(i)));
total_charge(:, i) = base_electricity_charge(:, i) + base_heat_charge(:, i);
figure(1);
plot(t_plot, power_demand(:, d_index)/1e3, t_plot, heat_demand(:, d_index)/1e3, 'r'); % kWh
legend('Electricity', 'Heat');
xlabel('Time (h)');
ylabel('Demanded power(kW)');
title('Electricity/Heat energy demand')
figure(2);
plot(t_plot, elec_tariff(:, d_index));
xlabel('Time (h)');
ylabel('Rate ($/kWh)');
title('Applicable electric rate');
%%
g = digraph(state_from, state_to, decided_costs(:, i));
[path, path_length] = shortestpath(g, 1, max(state_to), 'Method', 'acyclic');
[power_MGT(:, i), heat_MGT(:, i), mdot_MGT(:, i)] = extractPath(path, power_map, heat_map, fuel_map, SV_states);
path_cost(:, i) = path_length;
h1 = figure(3);
plot(t_plot, power_demand(:, d_index)/1e3, t_plot, power_MGT(:, i)/1e3, t_plot, (power_demand(:, d_index) - power_MGT(:, i))/1e3, 'LineWidth', 2); %in kW
legend('Power demand', 'CHP commitement', 'Utility commitement');
xlim([0, 25]);
ylim([0, Inf]);
grid on;
xlabel('Time [h]', 'FontSize', FS);
ylabel('Demand and commitment [kW]', 'FontSize', FS);
hold off;
h2 = figure(4);
plot(t_plot, heat_demand(:, d_index)/1e3, t_plot, heat_MGT(:, i)/1e3, t_plot, (heat_demand(:, d_index) - heat_MGT(:, i))/1e3, 'LineWidth', 2); %in Kw
legend('Heat demand', 'CHP commitement', 'Utility commitement');
xlim([0, 25]);
ylim([0, Inf]);
grid on;
xlabel('Time [h]', 'FontSize', FS);
ylabel('Demand and commitment [kW]', 'FontSize', FS);
%{
if save_fig
figname=['FC', num2str(fuel_index(i)), 'B', num2str(tariff_map(d_index,1)), 'D', num2str(tariff_map(d_index,2))];
print('-f1',[savepath, 'PH_Demand', figname],'-depsc');
print('-f2',[savepath, 'Tariff', figname],'-depsc' );
print('-f3',[savepath, 'Power_MGT', figname],'-depsc');
print('-f4',[savepath, 'Heat_MGT', figname],'-depsc');
end
%}
%Save just case study .fig files and eps
% if save_fig && fuel_index(i)==1
% figname=['FC', num2str(fuel_index(i)), 'B', num2str(tariff_map(d_index,1)), 'D', num2str(tariff_map(d_index,2))];
% savefig(h1,[savepath, 'Power_MGT', figname, '.fig']);
% savefig(h2,[savepath, 'Heat_MGT', figname, '.fig']);
% set(gcf, 'PaperPositionMode', 'auto');
% print('-f3',[savepath,'eps_noleg\', 'Power_MGT', figname],'-depsc')
% print('-f4',[savepath,'eps_noleg\', 'Heat_MGT', figname],'-depsc')
% end
new_demand(:, i) = subplus(power_demand(:, d_index)-power_MGT(:, i));
bought_elec(i) = sum(subplus(power_demand(:, d_index)-power_MGT(:, i)).*dt*joule2kWh.*elec_tariff(:, d_index)); %in $
sold_energy(i) = sum(subplus(-1.*(power_demand(:, d_index) - power_MGT(:, i))).*dt*joule2kWh.*elec_tariff(:, d_index)); %in $
bought_fuel(i) = sum(mdot_MGT(:, i)*dt*price_kg_f(fuel_index(i))); %in $
bought_heat(i) = sum(subplus(heat_demand(:, d_index)-heat_MGT(:, i)).*dt*joule2kWh.*heat_tariff(fuel_index(i))); %in $
MGT_cost(i) = bought_elec(i) - sold_energy(i) + bought_fuel(i) + bought_heat(i); %in $
[MGT_PDC(i), MGT_IDC(i), ut_PDC(i), ut_IDC(i), FC(i)] = generateDemandCharges(d_index, dt, power_demand(:, d_index), new_demand(:, i));
% Absolute savings:
% savings(:,i)=total_charge(:,i)-path_cost(:,i);
savings(i) = total_charge(i) - MGT_cost(i);
1;
end
disp(['All data saved to folder ', savepath]);
%% Save economic data
savepath = '..\Data';
save([savepath, 'econ_data.mat'], 'elec_tariff', 'heat_tariff', ...
'power_demand', 'heat_demand', 'power_MGT', 'heat_MGT', 'MGT_cost', ...
'electric_energy_kWh', 'heat_energy_kWh', 'base_electricity_charge', ...
'base_heat_charge', 'total_charge', 'bought_elec', 'sold_energy', ...
'bought_heat', 'bought_fuel', 'MGT_cost', 'savings', 'new_demand', ...
'FC', 'ut_PDC', 'ut_IDC', 'MGT_IDC', 'MGT_PDC', 'mdot_MGT');