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

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 = '[e-mail]'; % 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)

Answer 1

[Jasvin]

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!