How to put the desired value for LPSP in MATLAB code for optimizing a PV-battery system using PSO algorithm?

19 次查看（过去 30 天）
Trainee Engineer
Trainee Engineer 2023-1-6
回答： Harsh Sharma 2024-9-22
Ran in:
I want to optimize a PV-battery system using particle swarm optimization (PSO) algorithm for reducing the net present cost (NPC) and increasing reliability (LPSP equal to almost 0). I have no idea how to put the desired value for LPSP, for e.g., 1%, 5%, 10%, such that the program converges to that LPSP value as well as minimizing the NPC at the same time. The MATLAB code for the main PSO and the objective function are provided below.
Main_PSO.m
clc
clear
close all
%% Problem
nVar = 3; % number of variables
VarMin = [0 0 0]; % lower bound of variable
VarMax = [1000 1000 1000]; % upper bound of variable
%% PSO parameters
MaxIter = 100; % max number of iterations
nPop = 50; % population size
w = 1; % inertia
d = 0.99; % damping ratio of the inertia
c1 = 2; % acceleration 1
c2 = 2; % acceleration 2
%% Initial
x0.position = [];
x0.velocity = [];
x0.fitness = [];
x0.best.position = [];
x0.best.fitness = [];
x = repmat(x0, nPop, 1); % make a population
global_best.fitness = inf;
% generate initial population
for i = 1:nPop
% generate random solutions
for k = 1:nVar
x(i).position(k) = unifrnd(VarMin(k), VarMax(k));
end
x(i).velocity = zeros([1 nVar]); % initial velocity
x(i).fitness = objective_function(x(i).position); % calculate the fitness
x(i).best.position = x(i).position; % update the local best
x(i).best.fitness = x(i).fitness; % update the local best
if x(i).best.fitness < global_best.fitness
global_best = x(i).best;
end
end
Unrecognized function or variable 'objective_function'.
B = zeros(MaxIter, 1); % save the best fitness in each iteration
%% Main program
for j = 1:MaxIter
for i = 1:nPop
x(i).velocity = w*x(i).velocity + c1*rand([1 nVar]).*(x(i).best.position - x(i).position)...
+ c2*rand([1 nVar]).*(global_best.position - x(i).position); % update velocity
x(i).position = x(i).position + x(i).velocity; % update position
% check the range
for k = 1:nVar
if x(i).position(k) < VarMin(k)
x(i).position(k) = VarMin(k);
end
if x(i).position(k) > VarMax(k)
x(i).position(k) = VarMax(k);
end
end
x(i).fitness = objective_function(x(i).position);
if x(i).fitness < x(i).best.fitness
x(i).best.position = x(i).position; % update the personal best
x(i).best.fitness = x(i).fitness;
if x(i).best.fitness < global_best.fitness
global_best = x(i).best; % update the global best
end
end
end
w=w*d; % update the damping ratio
% save best fitness
B(j)=global_best.fitness;
disp(['Iteration ' num2str(j) '; Best fitness = ' num2str(B(j))]);
plot(B(1:j, 1), 'r.'); drawnow
end
objective_function.m
function NPC = objective_function(x)
%% variables
Npv = x(1); % number of PV panels
Nbat = x(2); % number of batteries
Ninv = x(3); % number of inverters
%% economics
lf = 25; % lifetime of PV system
I = 12.1/100; % inflation rate
d = 8/100; % discount rate
k = (1+I)/(1+d); % constant
%% PV specifications
Pp = 300; % rated power of PV module (W)
G_STC = 1000; % solar radiation incident on PV at STC (W/m2)
Tcoeff = -0.41; % temperature coefficient of power (%/oC)
Tc_STC = 25; % PV cell temperature under STC (oC)
fpv = 80/100; % PV derating factor
Tc_NOCT = 45; % nominal operating cell temperature (oC)
Ta_NOCT = 20; % ambient temperature at which the NOCT is defined (oC)
G_NOCT = 800; % solar irradiation at which the NOCT is defined (W/m2)
n_mp_STC = 18.3/100; % MPP efficiency under STC
const = 0.9; % product of solar transmittance and absorptance (Duffie and Beckman, 1991)
%% battery specifications
Qbat = 1000; % battery capacity (Wh)
sd = 0.0125/100; % battery self-discharge rate per hour
n_bat = 95/100; % charge efficiency of battery
Pbat_min = 0; % minimum battery capacity (Wh)
Pbat_max = Nbat*Qbat; % maximum battery capacity (Wh)
%% converter specifications
n_conv = 80/100; % converter efficiency
n_inv = 80/100; % inverter efficiency
%% weather data
A = readtable('weatherfile.xlsx'); % read weather file
Ta = A.(2); % ambient temperature (oC)
G = A.(7); % global solar irradiance (W/m2)
%% load profile
B = readtable('loaddemand.xlsx'); % read load demand file
L = B.(2); % hourly load demand (W)
Pload = repmat(L, 365, 1); % hourly load demand for one year (W)
%figure
%plot(Pload);
for t = 1:8760
Pload(t) = Pload(t);
end
Pload_total = sum(Pload);
%% PV modelling
Tc = (Ta+((Tc_NOCT-Ta_NOCT).*(G./G_NOCT).*(1-((n_mp_STC*(1-((Tcoeff/100)*Tc_STC)))/const))))./(1+((Tc_NOCT-Ta_NOCT).*(G./G_NOCT).*(((Tcoeff/100)*n_mp_STC)/const))); % PV cell temperature in the current time step (oC)
Ppv = (Npv*Pp).*(G./G_STC).*(1+((Tcoeff/100).*(Tc-Tc_STC))).*fpv; % output power of PV array in the current time step (W)
Ppv_total = sum(Ppv); % annual output power of PV array (W)
%% battery modelling
Pbat_state = Nbat*Qbat; % available energy in the storage at the beginning of the time step (Wh)
for t = 1:8760
if Ppv(t) > Pload(t)
Pbat(t) = (Pbat_state*(1-sd))+(((Ppv(t)*n_conv)-(Pload(t)/n_inv))*n_bat);
if Pbat(t) > Pbat_max
excess = Pbat(t)-Pbat_max;
deficiency = 0;
Pbat(t) = Pbat_max;
Pbat_state = Pbat_max;
else
excess = 0;
deficiency = 0;
Pbat_state = Pbat(t);
end
elseif Ppv(t) < Pload(t)
Pbat(t) = (Pbat_state*(1-sd))-(((Pload(t)/n_inv)-(Ppv(t)*n_conv))/n_bat);
if Pbat(t) < Pbat_min
excess = 0;
deficiency = Pbat_state-Pbat(t);
Pbat(t) = Pbat_min;
Pbat_state = Pbat_min;
else
excess = 0;
deficiency = 0;
Pbat_state = Pbat(t);
end
end
LPS(t) = Pload(t)-(Ppv(t)+Pbat(t));
end
Pbat = Pbat';
%% capital expenditures
Cpv = 2000; % capital cost of PV ($/kW)
Cbat = 376; % capital cost of battery ($/kWh)
Cinv = 700; % capital cost of inverter($/kW)
CAPEX = (Npv*(Pp/1000)*Cpv)+(Nbat*Cbat)+(Ninv*Cinv);
%% operation & maintenance expenditures
CM_pv = 33; % annual maintenance cost of PV ($/kW)
CM_bat = 5.17; % annual maintenance cost of battery ($/kW)
CM_inv = 0; % annual maintenance cost of inverter ($/kW)
OPEX = (Npv*(Pp/1000)*CM_pv)+(Nbat*CM_bat)+(Ninv*CM_inv); % total OPEX ($)
PW_OPEX = OPEX*k*((1-(k.^lf))/(1-k)); % present worth of OPEX ($)
%% replacement cost
% battery
for n = 5:5:15
PW_bat(n) = Nbat*Cbat*(k.^n);
end
PW_bat = sum(PW_bat);
% inverter
PW_inv = Ninv*Cinv*(k.^15);
% present worth of replacement cost
PW_R = PW_bat+PW_inv;
%% penalty factor
PF = 5.6; % in $/kWh
%% reliability factor
LPS_total = sum(LPS); % loss of power supply
LPSP = LPS_total/Pload_total; % loss of power supply probability
%% net present cost
NPC = CAPEX+PW_OPEX+PW_R+(PF*(LPS_total/1000));
All functions in a script must be closed with an 'end'.

请先登录，再回答此问题。

回答（1 个）

Harsh Sharma
Harsh Sharma 2024-9-22
Hi Trainee Engineer,
As per my understanding of the question, you want the program to converge to a target LPSP value while using the PSO algorithm to reduce NPC.
You can achieve this by adding a penalty to the NPC value in your “objective_function” based on the difference in your target LPSP value and the actual LPSP value. Below is the code to achieve the same, which you can add at the end of the “objective_function”.
LPS_total = sum(LPS);
LPSP = LPS_total / Pload_total;
% LPSP target as 1%
LPSP_target = 0.01;
% Calculate penalty
Penalty = 10000 * abs(LPSP - LPSP_target);
% net present cost with penalty
NPC = CAPEX + PW_OPEX + PW_R + (PF * (LPS_total / 1000)) + Penalty;
I hope this helps! May your next bug be an easy fix.

类别

在 Help Center 和 File Exchange 中查找有关 Programming 的更多信息

标签

产品


版本

R2016b

Community Treasure Hunt

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

Start Hunting!

Translated by