Optimisation of energy system

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

[derly]

Hi Hector. Did you finish your thesis? Where can I check the results?

Answer 1

[Alan Weiss]

Perhaps you will find some inspiration in the example Optimal Dispatch of Power Generators: Problem-Based.

Alan Weiss

MATLAB mathematical toolbox documentation