-
Notifications
You must be signed in to change notification settings - Fork 1
/
MC_HEATPOWER.m
162 lines (150 loc) · 7.39 KB
/
MC_HEATPOWER.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
154
155
156
157
158
159
160
161
162
function [RES] = MC_HEATPOWER(X, D)
% The Generators Allocation Matrix
% the Data structure (contains 2 structures 1 for the heat and one for the
% electrci network)
Wea=D.Weather; % the electric network data
Del=D.Del; % the electric network data
Dth=D.Dth; % the heat network data
% The decision variables
x_el=X.x_el; % the DGs atrucutre (PV,WT,EV,ST) in the electric network
x_th=X.x_th; % the DGs in the heat network
%% Define a Distributed Generator Object
% DGmdl=DG_module('D',Del); # --> not yet available, need flexible constructor
DGmdl=DG_module();
%% Scenarios cycle
[pdfENS, pdfCg,NETPav,NETuP, DGPu, ef] = deal(zeros(Del.SCn,1));
[nENS] = deal(zeros(Del.SCn,Del.NDn));
tic
% random initial time state
mnth = randi(12); % the random month of the year
day = 1;
td = randi(24); % random hour of the day td
%% Generate weather History
% ITot clear sky irradiance
% ws is the wind speed profile [m/s]
% Beta distributed solar irradiance [kW/m^2] the total 'real' irradiance
long=-3.27; % Barry Island UK Longitude
lat=51.3; % Barry Island UK Latitude
[ITot,ws,~,~]=Weather_Simulator(Del.SCn,long,lat,mnth*day);
s=betarnd(Wea.PVa(1) , Wea.PVb(1),size(ITot)).*ITot;
Weather.s=s; % save weather samples
Weather.ws=ws;% save weather samples
for i = 1:Del.SCn
% try a time-sequential simulation of the system state
td=td+1; % add 1 h
if td>24; td=1;day=day+1; end
if day>30; day=1; mnth=mnth+1; end
if mnth>12; mnth=1; end
% mnth = randi(12); % the random month of the year
% td = randi(24); % random hour of the day td
%% The Meccanical State of the component
% to include failures change a state from 1 to 0
FDmecst=ones(Del.FDn,1); % meccanical state of the cables (1= healthy)
NETmecst=ones(Del.NDn,Del.DGtn+1); % meccanical state of the generators (1= healthy)
%% A simple failure model
D.GenFailRates=[0.01 0.05 0.05 0.05 0.05]; % failure probability generators (probability of 1 --> 0 in 1-h)
D.FDFailRates= 0.01; % failure probability links
D.GenRecoveryRates= [0.6 0.3 0.3 0.3 0.3]; % recovery probability generators (probability of 0 --> 0 in 1-h)
D.FDRecoveryRates= 0.8; % recovery probability links
[FDmecst,NETmecst]=MarkovFailure(FDmecst,NETmecst,D);
MecStates.FDmecst(i,:)=FDmecst;
MecStates.NETmecst(i,:)=NETmecst(:);
%% Random electric and thermal loads
% electric loads
LDel = normrnd(Del.muLDel(:,td).*1.5, Del.stdLDel(:,td)); % in [KWel]
LDel(LDel<0)=0;
% thermal loads
LDth = normrnd(Dth.LDth(:,td).*1.5,(Dth.LDth(:,td).*1.5)/20); % expressed in [KWth]
LDth(LDth<0)=0;
LDtots.el(i)=sum(LDel);
LDtots.th(i)=sum(LDth);
%% Main supply generators
MSpow=Del.MSmu; % this can be the control variable
% MSpow = normrnd(Del.MSmu, Del.MSmu/10); % random power generated by the main sources
MSPav = MSpow.*NETmecst(:,1).*x_el(:,1);
% random external temperature conditioanl to month
Text=normrnd(Wea.muText(mnth),Wea.stdText(mnth)); % external environmental conditions
Weather.Text(i)=Text;
%% Power DGs
% PV
PVPav=DGmdl.Power_PV(s(i),Text,Del.NDn); % available power from PV pannles
PVPav = PVPav.*NETmecst(:,2).*x_el(:,2);
% WT
WPav=DGmdl.Power_WT(ws(i),Del.NDn); % available power from wind turbines
WPav = WPav.*NETmecst(:,3).*x_el(:,3);
% ST
STPav=DGmdl.Power_ST(Del.NDn); % available power from storages
STPav = STPav.*NETmecst(:,5).*x_el(:,4);
%% Heat generated by Heat Pumps (function of the extenral tamperature)
Twater=35; % target water temperature
% [P_onoffhp,COPdc_onoffhp]=OnOff_HP(Text,Twater);
[P_onoffhp,COPdc_onoffhp]=DGmdl.Power_HP_OnOff(Text,Twater);
% [P_idhp,COPdc_idhp]=ID_HP(Text,Twater,psi);
% Ncomp=randi(2); [P_mchp,COPdc_mchp]=MC_HP(Text,Twater,Ncomp);
HPav=P_onoffhp.*x_th;
Lel_onoffhp=HPav./COPdc_onoffhp;
Lel_onoffhp(isnan(Lel_onoffhp))=0;
LDel=LDel+Dth.LumpedThLoadsMatrix*Lel_onoffhp;
%% Simple Power-Heat-Balance
% solve flow-pression equations for the heat network
NeededPowTh=(LDth-HPav); % the thermal power needed is the load - HP available from the Heat Pumps
NeededPowTh(NeededPowTh<0)=0; % if there is a thermal production excedence, then it goes wasted!!
% assume 100% efficent (if is an electric backup system is quite colose to 100%)
LDelbackup= (Dth.LumpedThLoadsMatrix*(NeededPowTh));
LDel=LDel+ LDelbackup./Del.MScm; % MScm conversion efficiency
LDtots.th_to_elbackup(i)=sum(LDelbackup./Del.MScm);
LDtots.LumpedThLoads(i)=sum(Dth.LumpedThLoadsMatrix*Lel_onoffhp);
%% Run an OPTIMAL Power Flow, Here we have the Network Solver
NETPavM = [MSPav PVPav WPav STPav].*[ones(Del.NDn,1) Del.pckg(ones(Del.NDn,1),[1:2,4])]; %Matrix of electric power Available
NETPav(i) = sum(sum(NETPavM)); % Total available power produced
DGPp(i,:) = sum(NETPavM(:, 2:end));% Total production for each generator type
[Pgenerators(i,:) ,pdfENS(i), NETminCo, NETuP(i), DGPu(i), nENS(i,:), ef(i)] = OPF(x_el, NETPavM, LDel, FDmecst, Del);
%% Cost Model
%ep= electricity price is a function of the electric load
%NETCi= investment costs (fixed for x) prorated in thn years
%NETminCo= operational cost of the electrical network
%inc_th= incentives for producing clean thermal power (ASHP 2.61p/kWh incentive in pounds 2017 UK)
%Del.Inc= incentives for producing clean electric power in p/kWh
inc_th=0.0261;
ep = Del.Eprice(sum(LDel));
DGCi = Del.DGCi.*Del.pckg;
NETCi = (sum(sum(x_el).*[Del.MSCi DGCi([1:2,4])])+sum(x_th)*Dth.Cost_onoff)/(365*24*Del.thn);
HPPowused=sum(LDth-NeededPowTh);
NetGainDG_el=(Del.Inc + ep)*(DGPu(i))*Del.ts; % gain from producing renwable electric Energy
NetGainDG_th= inc_th*HPPowused*Del.ts; % gain from producing renwable electric Energy
pdfCg(i) = NETCi + NETminCo + - NetGainDG_el - NetGainDG_th;
% try to add a cost of the energy not supplied (which must be bought from the grid at a cost ep)
pdfCg(i) = pdfCg(i) + pdfENS(i)*ep;
%% Incentives 2017 UK
% ASHP 2.61 p/kWh
% NETPpM = [MSPp PVPp WPp EVPp STPp].*[ones(D.NDn,D.MSn) D.pckg(ones(D.NDn,1),:)];
% NETPp(i) = sum(sum(NETPpM));
% DGPavM = NETPavM(:, 2:end);
% DGPav(i) = sum(sum(DGPavM));
% TLD(i) = sum(LDel);
end
pdfCg(isnan(pdfCg))=[];% get rid of the NaN simulations
pdfENS(isnan(pdfENS))=[];% get rid of the NaN simulations
Pgenerators(isnan(pdfENS),:)=[];
RES.pdfPgenerators=Pgenerators; % power
RES.NETPav=NETPav;
RES.pdfnENS=nENS; % Enrgy not supplied per node
RES.pdfDGP_RES=DGPu; % Total renewable power electric used
RES.pdfNETuP=NETuP; % real used power in the net
RES.pdfDGPp=DGPp; % real used power of the dgs
RES.pdfLDel=LDtots; % total electric thermal and combined loads
RES.pdfENS=pdfENS;
RES.pdfCg=pdfCg;
RES.Weather=Weather; % samples of the weather
RES.MecStates=MecStates;% save mechanical states
% Some statistics
RES.EENS=mean(pdfENS);
RES.ECg=mean(pdfCg);
RES.CovENS=std(pdfENS)/RES.EENS;
RES.CovCg=std(pdfCg)/RES.ECg;
Alpha=0.95*100;
RES.p95ENS=prctile(pdfENS,Alpha);
RES.p95Cg=prctile(pdfCg,Alpha);
RES.Supequantile95_ENS=mean(pdfENS(pdfENS>=prctile(pdfENS,Alpha)));
RES.Supequantile95_ENSCg=mean(pdfCg(pdfCg>=prctile(pdfCg,Alpha)));
end