Optimisation of energy system

16 views (last 30 days)
Hector Aguirre
Hector Aguirre on 26 Jul 2021
Commented: derly on 25 Jun 2023
Hello,
For my Masters thesis I have developed the following script that finds the optimal technology mix and size that minimise the OPEX of the energy system of a hotel while meeting the electricity and heating demand at all times (all half-hour periods in a year).
The problem is that I now want to add a battery as another variable (and optimise its capacity) to allow the system to store in it any excess generation at any time to discharge it when generation is lower than demand, while conserving the objective function and constraints. Any idea of how I can do this?
Thanks in advance,
Hector
clear
% Load Demand Profile
load('DemandProfile.mat')
load('WeatherDataFayoum.mat','GTilt','vwind','Tamb')
% PARAMETERS
% PV
nompv = 0.3; % kW
areapv = 1.6; % m2
co2pv = 6e-3; % g CO2/Wh
effpv = 0.19; % elect efficiency
% PVT
nompvt = 0.2; % kW
areapvt = 1.2; % m2
co2pvt = 17e-3; % g CO2/Wh
effpvt = 0.15; % elect efficiency
effpvth = 0.48; % thermal efficiency
% ST
nomst = 1.4; % kW
areast = 2.14; % m2
co2st = 17e-3; % g CO2/Wh
effst = 0.64; % thermal efficiency
% WT
nomwt = 2.4; % kW
areawt = 10.87; % m2
rhoair = 1.225; % kg/m3
co2wt = 4e-3; % g CO2/Wh
effwt = 0.35;
% CHP
co2CHP = 32e-3; % g CO2/Whe
effCHP = 0.26; % elect efficiency
effCHPh = 0.62; % thermal efficiency
% Biomass
massbio = 0.2036; % kg/kWh
nombio = 15; % kW
co2bio = 32; % g CO2/kWh
effbio = 0.93; % thermal efficiency
for i=1:17520
% Demand data for different resolutions
demandelec(i) = (Elec(i)+Cooling(i))/0.5; % W
demandheat(i) = Heating(i); % W
% Calculate weather conditions at different resolutions
irradiance(i) = GTilt(i);
wind(i) = vwind(i);
temp(i) = Tamb(i);
% CONSTRAINTS: Meet electricity and heating demand at all times
A(i,:) = -[irradiance(i)*effpv*0.5,irradiance(i)*effpvt*0.5,0,0.5*rhoair*0.5*effwt*areawt*(wind(i)).^3,effCHP*0.5,0];
A(i+17520,:) = -[0,irradiance(i)*effpvth*0.5,irradiance(i)*effst*0.5,0,effCHPh*0.5,0.5*nombio*1000*effbio];
b(i,:) = -demandelec(i);
b(i+17520,:) = -demandheat(i);
end
Aeq = [];
beq = [];
% VARIABLES
x0 = [0,0,0,0,0,0];
lb = [0,0,0,0,0,0];
ub = [inf,inf,inf,inf,inf,inf]; % To be determined from area limitations, etc.
% Function
[x,V] = fmincon(@fObjFunc,x0,A,b,Aeq,beq,lb,ub);
PV = x(1);
PVT = x(2);
ST = x(3);
WT = x(4);
CHP = x(5);
BIO = x(6);
OPEX = V;
% OPTIMISE
function f = fObjFunc(x)
% Unpack the input vector
x1 = x(1); x2 = x(2); x3 = x(3); x4 = x(4); x5 = x(5); x6 = x(6);
% Reenter Parameters
% PV
opexpv_fixed = 3.96; % pounds/m2/year (2% of capex )
opexpv_var = 0; % pounds/kwh
effpv = 0.19; % elect efficiency
% PVT
opexpvt_fixed = 9.58; % pounds/m2/year (1.8% of capex )
opexpvt_var = 0.00125; % pounds/kwh
effpvt = 0.15; % elect efficiency
effpvth = 0.48; % thermal efficiency
% ST
opexst_fixed = 12.74; % pounds/m2/year (1.8% of capex )
opexst_var = 0.00125; % pounds/kwh (1250 L/MWh , 0.1 pence/L )
effst = 0.64; % thermal efficiency
% WT
opexwt_fixed = 138.84; % pounds/unit/year (5.2% of capex )
opexwt_var = 0.002; % pounds/kWh
areawt = 10.87; % m2
rhoair = 1.225; % kg/m3
effwt = 0.35;
% CHP
opexCHP_fixed = 0.073; % pounds/W/year
opexCHP_var = 0.005; % pounds/kWh
effCHP = 0.26; % elect efficiency
effCHPh = 0.62; % thermal efficiency
% Biomass
opexbio_fixed = 243.75; % pounds/unit/year (5% of capex )
opexbio_var = 0.003; % pounds/kWh
effbio = 0.93; % thermal efficiency
nombio = 15; % kW
% INDICATOR CALCULATIONS
load('WeatherDataFayoum.mat','GTilt','vwind','Tamb')
% Calculate yearly generation
for i=1:17520
% Calculate weather conditions at different resolutions
irradiance(i) = GTilt(i);
wind(i) = vwind(i);
% Generation (kWh)
GenPV(i) = x1 * irradiance(i) * effpv * 0.5/1000;
GenPVTe(i) = x2 * irradiance(i) * effpvt * 0.5/1000;
GenPVTh(i) = x2 * irradiance(i) *effpvth * 0.5/1000;
GenST(i) = x3 * irradiance(i) * effst * 0.5/1000;
GenWT(i) = (x4) * 0.5 * rhoair * effwt * areawt * (wind(i).^3) * 0.5/1000;
GenCHPe(i) = x5 * effCHP * 0.5/1000;
GenCHPh(i) = x5 * effCHPh * 0.5/1000;
GenBIO(i) = x6 * nombio * effbio * 0.5;
% Calculate OPEX
OpexPV = opexpv_fixed * x1 + opexpv_var * sum(GenPV);
OpexPVTe(i) = opexpvt_fixed * x2 + opexpvt_var * sum(GenPVTe);
OpexST(i) = opexst_fixed * x3 + opexst_var * sum(GenST);
OpexWT(i) = opexwt_fixed * x4 + opexwt_var * sum(GenWT);
OpexCHPe(i) = opexCHP_fixed * x5 + opexCHP_var * sum(GenCHPe);
OpexBio(i) = opexbio_fixed * x6 + opexbio_var * sum(GenBIO);
end
OPEXe = sum(OpexPV) + sum(OpexPVTe) + sum(OpexWT) + sum(OpexCHPe);
OPEXt = sum(OpexST) + sum(OpexBio);
% Objective function
f = OPEXe + OPEXt;
end
  1 Comment
derly
derly on 25 Jun 2023
Hi Hector. Did you finish your thesis? Where can I check the results?

Sign in to comment.

Answers (1)

Alan Weiss
Alan Weiss on 27 Jul 2021
Perhaps you will find some inspiration in the example Optimal Dispatch of Power Generators: Problem-Based.
Alan Weiss
MATLAB mathematical toolbox documentation

Categories

Find more on Thermal Analysis in Help Center and File Exchange

Products


Release

R2021a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!