How do I optimize two scalar variables for an output that relies on multiple different arrays?

3 次查看（过去 30 天）
Craig Atkinson
Craig Atkinson 2022-10-31
回答： Jasvin 2023-9-21
I have meterological data collected to model a photovoltaic system. I am trying to optimize tilt (Beta) and orientation/azimuth (Psi) of an arbitrary panel for any given location. The output being optimized is the median monthly energy production. I am doing this correctly right now with a loop, but because this loop relies on the tilt being an integer, and my azimuth can only be 180 or 360 degrees. It is simply choosing the orientation based on the latitude, and then calculating the potential energy production for each value of Beta. I have read the documentation for meshgrid(), ndgrid(), and ga() but I am not sure how any of these would work because all other vectors are raw data and must be in-time with each other. Is there a more elegant way to do this? This is not the first time I've run into a problem like this
I am only interested in the area between the two long comment lines, ignore the rest of the messy prototype code. I have included the full program for the sake of clarity, however it relies on an API key from the NSRDB. They are free and you can find one here https://developer.nrel.gov/signup/
clear
clc
%% INPUTS
geobasemap(geoaxes,'streets')
enableDefaultInteractivity(gca);
[lat,lon] = ginput(1);
geoscatter(lat,lon,'filled','b');
if exist('lat','var')
close gcf
end
lat = num2str(lat);
lon = num2str(lon);
if str2double(lon) < -30 % Used to index the better satellite data
index = 3;
else
index = 1;
end
In.api_key_NSRDB = '**********************************'; % Needs legitimate API key
In.wkt = append('POINT(',lon,'+',lat,')');
In.locate_base_NSRDB = 'https://developer.nrel.gov/api/solar/nsrdb_data_query.json?';
In.locate_url_NSRDB = [In.locate_base_NSRDB 'api_key=' In.api_key_NSRDB '&wkt=' In.wkt];
In.options = webread(In.locate_url_NSRDB);
In.reference_url = In.options.outputs(index).links(end).link;
In.interval = num2str(In.options.outputs(index).links(end).interval);
In.year = num2str(In.options.outputs(index).links(end).year);
In.utc = 'true';
In.email = 'paxejo9693@dicopto.com'; % random temp email
In.download_base_NSRDB = extractBefore(In.reference_url,"?");
realurl = [In.download_base_NSRDB '?' 'names=' In.year '&wkt=' In.wkt...
'&interval=' In.interval '&api_key=' In.api_key_NSRDB...
'&email=' In.email '&utc=' In.utc];
Out.SolarDataTable = webread(realurl,weboptions('Timeout',60,'ContentType','table'));
SolarData = table2struct(Out.SolarDataTable);
%% Assign each column of raw data to the corresponding header
Year = [SolarData.Year];
Month = [SolarData.Month];
Day = [SolarData.Day];
Hour = [SolarData.Hour];
Minute = [SolarData.Minute];
Temp = [SolarData.Temperature];
ClearskyDHI = [SolarData.ClearskyDHI];
ClearskyDNI = [SolarData.ClearskyDNI];
ClearskyGHI = [SolarData.ClearskyGHI];
DHI = [SolarData.DHI];
DNI = [SolarData.DNI];
GHI = [SolarData.GHI];
Zenith = [SolarData.SolarZenithAngle];
Albedo = [SolarData.SurfaceAlbedo];
WindSpeed = [SolarData.WindSpeed];
%% Create continuously increasing time vectors so that plots don't overlap
Period = str2double(In.interval)/60;
TimeFull = 0:Period:Period*(numel(Hour) - 1);
DOY = (TimeFull/24);
%% Analyze and plot the data
Phi = str2double(lat); % lattitude in degrees
Theta = str2double(lon); % longitude in degrees
% Scalar values
E0 = 1000; % reference irradiance
T0 = 25; % reference temperature celcius
U0 = 25; % constant heat transfer component
U1 = 6.84; % convective heat transfer component
Eta = 0.1707; % efficiency of solar panel
Absorption = 0.95; % absorption coefficient
PanelWattage = 335; % nameplate wattage of panel
Gamma = -0.004; % temperature degredation coefficient (-0.4%/°C)
E_sc = 1367; % solar constant
PanelSize = 1.96; % ZNshine panel area
PanelCount = 20; % Example grid size
DerateSystem = 0.86; % Derate taken from PVwatts v5
InverterEfficiency = 0.96; % Efficiency taken from PVwatts v5
Loss = DerateSystem*InverterEfficiency; % Total loss (really 1 - Loss)
% Vector values (1x8760)
B = 360/365 .* (DOY-81);
C = (360.*(DOY./365));
Ratio = (1.00011 + 0.034221.*cosd(C) + 0.00128.*sind(C)...
+ 0.000719.*cosd(2.*C) + 0.000077.*sind(2.*C)); % ratio between distance from sun and average distance from sun
EoT = 9.87.*sind(2.*B) - 7.53.*cosd(B) -1.5.*sind(B); % minutes
TC = 4.*Theta + EoT; % minutes
LT = Hour + Minute/60; % local time
LST = LT + TC/60; % local solar time
HRA = 15.*(LST - 12); % hour angle in degrees
Delta = 23.45.*sind(B); % declination angle in degrees
Alpha = 90 - Zenith; % elevation/altitude angle in degrees
Azi = real(acosd((sind(Delta).*cosd(Phi) - cosd(Delta).*sind(Phi).*cosd(HRA))./cosd(Alpha)));
Azimuth = 360 - Azi; % Azimuthal angle in degrees
idx = HRA > 0;
Azimuth(~idx) = Azi(~idx);
E_a = max(E_sc.*Ratio,0); % extraterrestrial irradiance
Ai = DNI./E_a; % anisotropy index
% This is where I am trying to optimize.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if Phi > 0
Psi = 180; % module orientation angle (180 is south)
else
Psi = 360;
end
Beta = 1; % Starting module tilt angle
Betamax = 90;
Betavec = zeros(Betamax,1);
while Beta <= Betamax
AOI(Beta,:) = acosd(cosd(Zenith).*cosd(Beta) + sind(Zenith).*sind(Beta).*cosd(Azimuth - Psi)); % angle of incidence
SVF(Beta,:) = (1 + cosd(Beta))./2; % sky view factor
E_Beam(Beta,:) = max(DNI.*cosd(AOI(Beta,:)),0); % Sandia beam irradiance model
% squared ratio of the mean distance between the earth and the sun and
% the actual distance between the earth and the sun
% Reindl sky diffuse irradiance
E_diffuse(Beta,:) = real(max(DHI.*(Ai.*cosd(AOI(Beta,:)) + (1-Ai).*(SVF(Beta,:))...
.*(1 + sqrt((DNI.*cosd(Zenith))./GHI).*(sind(Beta/2)).^3)),0));
E_albedo(Beta,:) = max(GHI.*Albedo.*(1 - SVF(Beta,:)),0); % Sandia albedo irradiance model
E_POA(Beta,:) = E_Beam(Beta,:) + E_diffuse(Beta,:) + E_albedo(Beta,:); % Sandia plane of array irradiance
E_Month(Beta,:) = trapz(reshape(E_POA(Beta,:),[],12));
OptimalParam(Beta,:) = median(E_Month(Beta,:));
Betavec(Beta) = Beta;
Beta = Beta + 1;
end
Opt.Psi = Psi; % Opt is just a structure to hold the final design parameters
Opt.Beta = Betavec(OptimalParam == max(OptimalParam));
Opt.AOI = acosd(cosd(Zenith).*cosd(Opt.Beta) + sind(Zenith).*sind(Opt.Beta).*cosd(Azimuth - Psi)); % angle of incidence
Opt.SVF = (1 + cosd(Opt.Beta))./2; % sky view factor
Opt.E_Beam = max(DNI.*cosd(Opt.AOI),0); % Sandia beam irradiance model
Opt.E_diffuse = real(max(DHI.*(Ai.*cosd(Opt.AOI) + (1-Ai).*(Opt.SVF)...
.*(1 + sqrt((DNI.*cosd(Zenith))./GHI).*(sind(Opt.Beta/2)).^3)),0)); % Sandia diffuse irradiance model
Opt.E_albedo = max(GHI.*Albedo.*(1 - Opt.SVF),0); % Sandia albedo irradiance model
Opt.E_POA = Opt.E_Beam + Opt.E_diffuse + Opt.E_albedo; % Sandia plane of array irradiance
Opt.E_Month = trapz(reshape(Opt.E_POA,[],12)); % 1x12 array
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Tm = Temp + Opt.E_POA./(U0 + U1.*WindSpeed); % David Faiman module temperature model
Tc = Temp + Opt.E_POA.*((Absorption.*(1-Eta))./(U0 + U1.*WindSpeed)); % PVSyst cell temperature model
P_dc = (Opt.E_POA./E0).*PanelWattage.*(1 + Gamma.*(Tc - T0)); % Output power (kW/m^2)
P_expected = P_dc.*Loss;
P_grid = P_expected*PanelCount;
figure(4)
subplot(2,1,1)
boxplot(transpose(E_Month),'PlotStyle','compact')
subplot(2,1,2)
bar(Opt.E_Month)

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

回答（1 个）

Jasvin
Jasvin 2023-9-21
Hi Craig,
If I understand your problem correctly, you have an optimization problem where you are trying to optimize the value of the median monthly energy production which you are storing in the variable “OptimalParam(Beta, :)" which itself is dependent on several quantities that you are defining in the “for” loop but they are all in turn dependent on “Beta”, “Psi” and other constants from the data that you are importing.
You have two optimization variables, “Psi” and “Beta”. The constraints are that “Psi” can take only two values 180 and 360 and “Beta” is an integer variable that can take values in the range of 1 to 90.
Modelling discrete constraints like the one on “Psi” is a bit complicated but can be done. The following MATLAB Answers page describes how one can accomplish this:
https://www.mathworks.com/matlabcentral/answers/334289-how-can-i-perform-optimization-over-a-discrete-set-of-possible-values-in-matlab
The MATLAB Answer primarily links to this documentation guide, refer to the “Add Discrete Non-Integer Variable Constraints” section for the specific requirements:
https://www.mathworks.com/help/gads/solving-a-mixed-integer-engineering-design-problem-using-the-genetic-algorithm.html
Please note that you will also have to refer to the following MATLAB example,
openExample('globaloptim/steppedCantileverExample')
And examine the following dependency scripts to fully understand the scope of the example and for providing discrete integer constraints:
% For understanding how the fitness/objective function is defined
edit cantileverVolume.m
% For understanding how the discrete integer variables are mapped
edit cantileverMapVariables.m
edit cantileverConstraintsWithDisc.m
edit cantileverConstraints.m
The guide will also help you in modelling the integer constraint of the “Beta” variable.
Hope this helps!

类别

在 Help Center 和 File Exchange 中查找有关 Fixed-Wing Aircraft Creation with Objects 的更多信息

标签

产品


版本

R2021b

Community Treasure Hunt

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

Start Hunting!

Translated by