Looking at the highest values for each time of the day we get an asymmetrical shape with lower maximum values in the morning in June while we get less than expected maximum values around noon in December. This pattern is caused by two trees to the east and south of my home which cast a shade on the roof during different times of the day. The tree to the south is not too close to the house and only when the Sun’s elevation is very low in winter does it’s shade reach the roof.

In the first post on my solar panel dataset we did some exploration and found a number of interesting patterns. One of these was the effect of temperature on the efficiency of the panels which we investigated further in the second post . Another interesting pattern found while exploring the data was that for a number of months the maximal power throughout the day did not follow the expected bell shaped curve. Let’s load the data and visualize the maximal values per time of the day for June and December so you get an idea of what I’m talking about.

While I know that the trees are behind this pattern I was wondering if we could derive their location and shape from the data. If we could somehow get a value for the power loss for each timestamp in the dataset we should be able to couple this data to a position of the Sun at that time using the solarpos function from the maptools package. Then we could create a panoramic picture where each pixel is a position of the Sun throughout the year as seen by the panels. The solar positions that consistently yield too little power should then match the trees. This sounds like a fun idea, but to get there we first need to get an idea of the power loss per timestamp. We’ll have to create a model that gives us the theoretical maximum output for each timestamp and then calculate the difference with the observed values.

These are the Sun’s positions throughout the year on each round hour (the 10 minute intervals are a bit too much to plot). The resulting shapes are called analemmas and they are really interesting. Back in ancient times they gave astronomers headaches since they cause small deviations in the time read from sundials throughout the year. The issue was eventually resolved when they came up with the equation of time . I recommend a dive down this Wikipedia rabbit hole but for now, let’s get back to our solar panel data!

Before we dive into modelling the maximal power outputs let’s look at the solar positions we will bind the power loss data to. We get all solar positions for a full year in 10 minute intervals and filter out the values where the Sun is below the horizon.

Modelling maximal power output

We want a model that produces nice, symmetrical, bell shaped curves that are bigger in summer and smaller in winter. We could go for an approach where we simply try to create the expected curves using, for example, a normal distribution and some seasonal variable. While this approach would certainly work, it would treat the underlying effects as a black box. Wouldn’t it be more interesting if we could actualy simulate the effects at play and let the expected curve emerge by combining them? To do this we have to think about what is behind the variation in power yield. In winter the Sun transfers less energy to the panels because of it’s lower elevation but what exactly causes this power loss?

Energy loss due to spreading Image you carry a flashlight in a dark room. When you point it directly at the wall (90° angle) you get a small circle that is intensely lit. However, if you stand close to the wall and point the flashlight somewhat sideways you light a bigger surface. The same energy emitted by the flashlight is now spread out and each square centimeter in the lit area now receives less. We can calculate this effect with cos(pi/2 - angle) . When we transform a 90° angle to radials we get pi/2 and when we enter this in the function we get the cosine of zero which is 1 or 100% spreading efficiency. In Belgium the Sun has a maximal elevation of 62° which has a spreading efficiency of 88%. During the winter solstice the Sun only reaches 15.3° at noon which has a spreading efficiency of only 26%. Fortunately, the panels themselves are somewhat inclined (20°) towards the south so that even in winter the angle of arrival is closer to 35° (57% efficiency). However, this inclination results in a more complex spreading effect because the panels now point in a particular direction. Both effects are estimated in the code below. # angles in degrees to angles in radials # the '_deg' part of the variable name shows it's in degrees deg_to_rad <- function(angle_deg){pi * angle_deg / 180} # vertical spreading effect get_vertical_spreading <- function(angle_deg){ cos(pi / 2 - deg_to_rad(angle_deg)) } # When the Sun points to the south we should add 20° to the # solar elevation to get the correct angle since the panels # are inclined by 20° towards the south. # However if the Sun were to shine from the north we should # detract 20° from the elevation. The negative cosine # should do just this for us. adjust_horizontal_angle <- function(angle_deg, azimuth_deg){ angle_deg + 20 * -cos(deg_to_rad(azimuth_deg)) }

Energy loss due to traveling through the atmosphere A second, less straightforward, effect that decreases the energy transfer from the Sun to the panels at lower solar elevations is the fact that the distance traveled through the atmosphere becomes larger. As the light rays make their way through this mass of air they often collide with the molecules and become scattered. It is this effect that makes sunsets so beautiful and makes the sky above above us appear blue instead of black. Some clever people have studied this effect in great detail and we can take their equations from Wikipedia to learn that at a 90° angle the energy per square meter is about 1051 Watts. At 62° (Belgium summer solstice) we get 1019 Watts/m² which is still quite good. However at 15.3° in winter we fall back to 628 Watts/m². While we lose a lot of energy due to this scattering effect a little bit may be recuperated. Some of the scattered light from rays passing overhead now reaches our panels. # calulate Air Mass (AM) according to Karsten & Young # https://en.wikipedia.org/wiki/Air_mass_(solar_energy)#math_I.1 # https://en.wikipedia.org/wiki/Solar_irradiance get_air_mass <- function(angle_deg){ 1 / (cos(deg_to_rad(angle_deg)) + 0.506 * (96.08 - angle_deg) ^ -1.6364) } # this function calculates how much power (W/m²) # is left in Sun rays reaching the surface # the only effect weakening solar power is air mass at this point # incoming solar power on top of atmosphere = 1366 W/m² get_solar_intensity <- function(angle_deg){ # extract from 90 to get correct z angle (see wikipedia) 1.1 * 1366 * 0.7 ^ get_air_mass(90 - angle_deg) ^ 0.678 } # Raw estimation of the scattered light reaching the panels get_scatterlight_intensity <- function(angle_deg){ # The solar intensity of rays passing overhead solar_intensity <- get_solar_intensity(angle_deg + 3) # The percentage of max energy (at 90°) that remains pct_intensity <- solar_intensity / get_solar_intensity(90) # Scattered light is a fraction of the lost energy scatter_intensity <- solar_intensity * (1 - pct_intensity) * 0.15 return(scatter_intensity) }

Putting it all together Now that we have created functions for the most important effects at play it’s time to put them together and create a master function that will estimate the maximal power yield of the solar installation when we feed it the position of the Sun. solar_position_to_max_power <- function(angle_deg, azimuth_angle_deg){ # how much of the solar intensity is left # after travelling through the air mass? solar_intensity <- get_solar_intensity(angle_deg) # light spreading out because of sideways spreading horizontal_angle_adjusted <- adjust_horizontal_angle(angle_deg, azimuth_angle_deg) # panels are inclined at 20° vertical_spreading <- get_vertical_spreading(horizontal_angle_adjusted) surface_power <- {solar_intensity * vertical_spreading + get_scatterlight_intensity(angle_deg)} # 7 panels with 1.6 m² panel_surface <- 1.6 * 7 # 16% efficiency of panels efficieny <- 0.18 # eventual power is the panel surface times # the power left in the Sun's rays on the panel times # the efficiency of the panels max_plant_power <- panel_surface * surface_power * efficieny return(max_plant_power) } We can now call this function on the dataset with solar positions for a full year which we created earlier on. solar_position_max_power_df <- solar_position_df%>% # calculate the maximal power of the solar panels per timestamp # we need this to check if the maximal observed power is much lower # than the maximal theoretical power mutate(max_power = solar_position_to_max_power(solar_elevation, solar_azimuth_angle))%>% filter(max_power >= 0) Now let’s bring both datasets together and calculate the ratio between the power observed and the estimated maximal power. combined_df <- solarplant_df%>% # add the solar plant data to the solar position data # standardise the time to 1970 to allow joing the solar position dataset mutate(timestamp_basic = lubridate::origin + (yday(timestamp) - 1) * 3600 * 24 + hour(timestamp) * 3600 + minute(timestamp) * 60)%>% right_join(solar_position_max_power_df, by = c("timestamp_basic" = "timestamp"))%>% replace_na(list(power = 0))%>% # efficiency is power / max power but can't be higher than 1 (100%) mutate(efficiency = ifelse(power / max_power > 1, 1, power / max_power), # some time variables we need later on hour = hour(timestamp) + minute(timestamp) / 60, week_date = round_date(timestamp_basic, unit = "week")) To check if our theoretical model delivers acceptable results let’s visualize both observed and estimated maximal power for some sunny days in winter and summer. combined_df%>% mutate(date = as_date(timestamp))%>% filter(date %in% as_date(c("2017-06-01", "2016-12-04")))%>% ggplot(aes(x = hour(timestamp) + minute(timestamp) / 60, y = max_power))+ geom_line()+ geom_point(aes(y = power))+ facet_wrap(~date)+ labs(y = 'Power (W)', x = 'Time of day (h)')+ theme_minimal()+ theme(text = element_text(colour = "grey40"), strip.text.x = element_text(colour = "grey40"), panel.grid.minor = element_blank(), panel.border = element_rect(colour = "grey40", size = rel(0.3), fill = NA)) As you can see our theoretical estimation fits the observed values quite well both in winter and summer. While not perfect, this estimation should now allow us to detect solar positions that never result in maximal efficiency values due to an obstacle blocking the light.