# Optimisation of energy system

9 views (last 30 days)
Hector Aguirre on 26 Jul 2021
Answered: Alan Weiss on 27 Jul 2021
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?
Hector
clear
% 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
wind(i) = vwind(i);
temp(i) = Tamb(i);
% CONSTRAINTS: Meet electricity and heating demand at all times
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
% Calculate yearly generation
for i=1:17520
% Calculate weather conditions at different resolutions
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

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

R2021a

### Community Treasure Hunt

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

Start Hunting!