Optimal Placement of Solar Panels#
About the project#
Duration: 20-50 hours (including typesetting final report in LaTeX)
Prerequisites: Multivariable calculus (vector fields, surface integrals), basic linear algebra, coordinate systems, basic Python programming
Python packages: NumPy, matplotlib, pandas, pvlib, scipy
Learning objectives: Model solar irradiance, calculate optimal panel angles, understand vector field flux through surfaces, apply numerical integration methods
Introduction#
A homeowner near DTU wants to install solar panels on his detached house with a flat roof. She thinks she should orient them towards the south but is unsure about the angle the panel should form with the roof. In this project, you are to find the optimal angle for the homeowner’s solar panels. The goal of the project is to develop a model and a Python program that, given the location on Earth (latitude and longitude), can calculate the optimal angle for the installation of solar panels under a set of simplifying assumptions. Data from the calendar year 2024 will be used for simulation. You will begin by modeling the movement of the sun in the sky using the Python package pvlib. Then, you will establish a model for the solar panel’s energy production, which does not account for clouds and the atmosphere. The sun is modeled as a vector field of parallel vectors. To calculate the energy production, you must be able to integrate the flux of the vector field through the solar panel surface. Initially, we assume the panel is flat, but you will subsequently generalize your model. One possibility for generalization is as follows: Find a suitable surface in the urban environment around DTU that can be equipped with thin, flexible solar panel films (like https://www.youtube.com/watch?v=TS9ADU0oc50) and calculate the energy production from such solutions.
Preparation#
From <a href=https://01002.compute.dtu.dk/_assets/textbook_mat1b_en.pdf#exe:max-direct-deriv-in-grad”>this book study:
Orthogonal projections, Section 2.3
Spherical coordinates (Example 6.6.1) and the horizontal coordinate system (https://en.wikipedia.org/wiki/Horizontal_coordinate_system)
Surfaces and normal vector, Section 7.1
Vector fields and flux, Sections 7.2 and 7.5
About the Project#
Project Purpose#
The overall goal of the project is to create a Python program that, given a specified location on Earth (latitude and longitude), can calculate the optimal angle for solar panels.
Optimality#
“Optimal” can mean many things in this context:
Greatest annual energy production. This is the optimality criterion we use unless otherwise specified.
Requirement for a certain daily minimum output (e.g., during winter).
Lowest total energy costs for the residence. This should include consumption patterns, selling energy to the grid, potential electric vehicle usage, etc.
Etc.
Theory: Model and Assumptions#
Solar energy is the conversion of energy from sunlight into electricity using solar cells (photovoltaic (PV) cells). The production of power by solar panels depends on the size of the installation (the number of panels), the individual efficiency of the panels, their placement, and, of course, weather conditions.
Solar irradiance, or solar irradiance, is the power per unit area (surface power density) received from the sun in the form of electromagnetic radiation in the instrument’s wavelength range. Solar irradiance is measured in watts per square meter (\(\mathrm{W/m^2}\)) in SI units. Solar radiation has different wavelengths, such as UV radiation, which collectively are called the spectrum of radiation.
The average annual solar radiation arriving at the top of the Earth’s atmosphere is about \(1361 \,\mathrm{W/m}^2\). When the sun’s rays have passed through the atmosphere, the total radiation (irradiance) in Denmark is approximately \(S_0 := 1100 \,\mathrm{W/m^2}\) on a clear day at noon. The spectrum also changes through the atmosphere; for example, the shortest UV radiation (below 280 nanometers) rarely reaches the Earth’s surface, but we will ignore such spectral changes.
With a good approximation, one can assume a linear relationship between the solar irradiance and the current output of the solar panel, and we will use this approximation. However, since the temperature in practice increases with increased irradiance, and the voltage decreases with increased temperature, the effect will not increase entirely linearly. The sun’s rays are modeled as a vector field \(\pmb{V}_S\) of parallel vectors of length \(S_0\). Let’s assume that the solar panel is described by a surface \(\mathcal{F} = \pmb{r}(\Gamma)\), where \(\pmb{r} : \Gamma \to \mathbb{R}^2\), \(\Gamma = [a_1,b_1] \times [a_2,b_2]\), is a parameterization of the surface. The energy production for a solar cell depends on the angle between the surface’s normal vector \(\pmb{n}_\mathcal{F}\) and the solar rays. If the solar rays are parallel to the normal vector, the solar power is maximum, while if the solar rays are perpendicular to the normal vector, the power is zero. The solar power per unit area at the point \(\pmb{r}(u,v)\), where \((u,v) \in \Gamma\), is the projection of the solar rays onto the surface’s unit normal vector:
The total effect is obtained by integrating this quantity over the entire surface, which corresponds to the flux of \(V\) through the surface \(\mathcal{F}\).
If \(\langle \pmb{V}, \pmb{n}_\mathcal{F}(u,v) \rangle < 0\), it means that the sun is shining on the backside of the solar panel, and thus the power should be set to zero (the power can never be negative).
Unless otherwise stated, we adopt the following standard assumptions:
The solar panel is flat and fixed in orientation.
The maximum solar irradiance is assumed to be \(S_0 = 1100 \,\mathrm{W/m^2}\). All vectors in the sun’s vector field \(\pmb{V} : \mathbb{R}^3 \to \mathbb{R}^3\) are parallel and have magnitude \(S_0\).
The power output of the solar panel depends linearly on the solar irradiance, modeled as the flux of the sun’s vector field \(\pmb{V}\) through the surface of the panel.
The efficiency \(\eta\) of the solar panel is defined based on the standard test condition (STC) of \(1000 \,\mathrm{W/m^2}\) irradiance. At this flux, the panel delivers
\[\begin{equation*} \frac{Wp}{L \cdot B} \,\mathrm{W/m^2}, \end{equation*}\]where \(Wp\) is the panel’s peak power (in watts), and \(L\) and \(B\) are its length and width (in meters), respectively.
Cloud cover and atmospheric disturbances are not modeled explicitly. Instead, we assume an average annual attenuation factor \(A_0 = 0.5\), which effectively halves the incoming solar energy:
\[\begin{equation*} S_0 A_0 = 550 \,\mathrm{W/m^2}. \end{equation*}\]Thus, we treat the sun’s vector field \(\pmb{V}\) as having constant direction and reduced magnitude \(A_0 S_0\).
In the final part of the project, there will be an opportunity to revisit these assumptions and investigate whether/how they can be replaced with better assumptions.
Initial Exercises#
Find through literature search the recommended angle at which solar panels are installed in Denmark. The angle is said to be zero degrees if the panel lies flat on the ground (or roof).
Select a type of solar panel. You can search for installers in Denmark and investigate which panels they typically use, or you can google “solar panel datasheet” or similar. The solar panel should be a standard panel (flat and not curved, for example). Find a datasheet for the chosen solar panel, and describe the panel’s size (using \(L\) for length and \(B\) for width) and indicate Wp/Pmax (referred to as max power, peak power watts, or similar) under the STC standard. Describe what the STC standard entails. Calculate \(Wp/(L B)\) (according to the list of standard assumptions above). Assume ideal conditions: the sun is perpendicular to the solar panel, and the solar irradiance is \(1100 W/m^2\) for an entire hour. How many joules and kilowatt-hours does the panel deliver in this hour? How many watt-hours is this per \(m^2\)?
Your final report should also include an introductory section where you describe solar energy and solar cells. You determine what exactly the section should contain, but the sources below can be used to find information. It is recommended to write this section first once you have progressed further with the project.
Since the panel is flat, the unit normal vector \(\pmb{u}_p\) is constant (and has length 1), where:
Provide a formula or expression for the flux through the surface expressed in terms of \(A_0\), \(\pmb{u}_p, \pmb{V}, L\), and \(B\). We define the flux to be zero if the angle between \(\pmb{u}_p\) and \(\pmb{V}\) is greater than \(\pi/2\) (90 degrees), as we do not want negative flux. Your expression should account for this.
Hint: Since the normal vector is constant, you can get rid of the integral sign.
The flux determines the instantaneous power for the panel for a specific direction of the sun’s rays \(\pmb{V}\) and for a specific panel orientation determined by \(\pmb{u}_p\). To find the total energy production for a solar panel, we need to be able to model the movement of the sun relative to the panel. This is modeled using a Solar Position Algorithm (SPA) from the Python package pvlib, but before we get to that point, we need to develop some tools in NumPy.
Provide SI units for \(\pmb{V}\), \(L\), \(B\), \(A_0\), flux, and energy. Specify the relationship between \(\mathrm{J}\) and \(\mathrm{kWh}\).
NumPy#
In the project, we need to be able to find minimum, maximum, and/or zero points as NumPy arrays, so we start by familiarizing ourselves with this. We consider a simple function \(f: [0, 2 \pi] \to \mathbb{R}\) given by \(f(x) = \cos(x)\) and want to determine its minimum, maximum, and zero points. However, in the project, we won’t have the actual function or its functional expression, but only a (long) array of function values, for example, \(f(t_n)\) for \(t_n = 2\pi n/N\) for \(n = 0,1, \dots, N-1\). Here we perceive \(\Delta t_n = 2 \pi/N\) as the distance between our time samples. In NumPy:
import numpy as np
import sympy as sym
import matplotlib.pyplot as plt
figuresize = (15, 7)
t = np.linspace(0, 2 * np.pi, 1000)
f = np.cos(t)
plt.figure(figsize=figuresize)
plt.plot(t, f, color='b', linestyle='-')
plt.xlabel('t')
plt.ylabel('cos(t)')
plt.show()
# print(t,f) # uncomment to see the values of t and f
We perceive t and f as vectors of variable and function values, respectively. NumPy has many built-in methods that can be useful, such as
f.max(), f.min(), f.argmax(), f.argmin(), t[f.argmax()], t[f.argmin()]
(1.0, -0.9999950553174459, 0, 499, 0.0, 3.138447916198812)
One can also inquire about which function values, for example, are less than \(-0.95\):
idx = f < -0.95
f[idx], t[idx]
# Or in one line
f[f < -0.95], t[f < -0.95]
(array([-0.95192731, -0.95383508, -0.95570513, -0.95753737, -0.95933173,
-0.96108814, -0.96280654, -0.96448685, -0.96612901, -0.96773295,
-0.96929861, -0.97082592, -0.97231483, -0.97376528, -0.97517722,
-0.97655057, -0.9778853 , -0.97918134, -0.98043865, -0.98165717,
-0.98283687, -0.98397768, -0.98507957, -0.9861425 , -0.98716641,
-0.98815128, -0.98909706, -0.99000371, -0.9908712 , -0.99169949,
-0.99248855, -0.99323836, -0.99394887, -0.99462007, -0.99525192,
-0.9958444 , -0.99639749, -0.99691116, -0.9973854 , -0.99782019,
-0.9982155 , -0.99857133, -0.99888765, -0.99916446, -0.99940175,
-0.99959951, -0.99975772, -0.99987639, -0.9999555 , -0.99999506,
-0.99999506, -0.9999555 , -0.99987639, -0.99975772, -0.99959951,
-0.99940175, -0.99916446, -0.99888765, -0.99857133, -0.9982155 ,
-0.99782019, -0.9973854 , -0.99691116, -0.99639749, -0.9958444 ,
-0.99525192, -0.99462007, -0.99394887, -0.99323836, -0.99248855,
-0.99169949, -0.9908712 , -0.99000371, -0.98909706, -0.98815128,
-0.98716641, -0.9861425 , -0.98507957, -0.98397768, -0.98283687,
-0.98165717, -0.98043865, -0.97918134, -0.9778853 , -0.97655057,
-0.97517722, -0.97376528, -0.97231483, -0.97082592, -0.96929861,
-0.96773295, -0.96612901, -0.96448685, -0.96280654, -0.96108814,
-0.95933173, -0.95753737, -0.95570513, -0.95383508, -0.95192731]),
array([2.83026365, 2.83655313, 2.8428426 , 2.84913208, 2.85542155,
2.86171103, 2.8680005 , 2.87428998, 2.88057945, 2.88686892,
2.8931584 , 2.89944787, 2.90573735, 2.91202682, 2.9183163 ,
2.92460577, 2.93089525, 2.93718472, 2.9434742 , 2.94976367,
2.95605315, 2.96234262, 2.9686321 , 2.97492157, 2.98121105,
2.98750052, 2.99379 , 3.00007947, 3.00636895, 3.01265842,
3.0189479 , 3.02523737, 3.03152684, 3.03781632, 3.04410579,
3.05039527, 3.05668474, 3.06297422, 3.06926369, 3.07555317,
3.08184264, 3.08813212, 3.09442159, 3.10071107, 3.10700054,
3.11329002, 3.11957949, 3.12586897, 3.13215844, 3.13844792,
3.14473739, 3.15102687, 3.15731634, 3.16360582, 3.16989529,
3.17618476, 3.18247424, 3.18876371, 3.19505319, 3.20134266,
3.20763214, 3.21392161, 3.22021109, 3.22650056, 3.23279004,
3.23907951, 3.24536899, 3.25165846, 3.25794794, 3.26423741,
3.27052689, 3.27681636, 3.28310584, 3.28939531, 3.29568479,
3.30197426, 3.30826374, 3.31455321, 3.32084268, 3.32713216,
3.33342163, 3.33971111, 3.34600058, 3.35229006, 3.35857953,
3.36486901, 3.37115848, 3.37744796, 3.38373743, 3.39002691,
3.39631638, 3.40260586, 3.40889533, 3.41518481, 3.42147428,
3.42776376, 3.43405323, 3.44034271, 3.44663218, 3.45292166]))
where we also specify the corresponding time values in the t vector.
Write Python code that finds all function values in
fin the interval \([-0.05, 0.05]\) and indicates the correspondingtvalues.
Write a Python function that can find both sign changes (zero crossings) of
f.
Check that the found zero crossings reasonably align with the exact zeros for \(f\). Your method should be simple and also work on “vectors” np.array([3, 0.5, 0.1,-0.5, -1, -0.5, 0.5, 2]) where we don’t know the function behind the function values. In particular, your method should not use any other knowledge about \(\cos\) (such as its differentiability) than the information you have in the f vector.
You can find inspiration in the bisection method (which you encountered in Mathematics 1a in the fall), or you can find inspiration in NumPy functions np.where, np.diff, np.sign, or the following Python plot:
plt.figure(figsize=figuresize)
plt.plot(t, np.abs(f), color="r") # linestyle='-', marker='.'
plt.xlabel("t")
plt.ylabel("|cos(t)|")
plt.show()
Solar Position Model#
Coordinate System#
To be able to calculate, for example, the movement of the sun and the projection of the sun’s rays onto the panel’s normal vector, we need to introduce (at least) an appropriate coordinate system. We choose to place a fixed center (origin) \((0,0,0)\) at the panel’s position. We then imagine the tangent plane to the Earth at the panel’s position. The normal vector to this tangent plane will be our \(z\)-axis. The \(x\)-axis and \(y\)-axis thus “span” the tangent plane. We define the coordinate system such that the \(x\)-axis points north and the \(y\)-axis points east. The attentive reader will have noticed that this is a left-handed coordinate system. However, this is the tradition, and the coordinate system is called the horizontal coordinate system (https://en.wikipedia.org/wiki/Horizontal_coordinate_system). So now we have a reference solar position system that focuses on the observer (the panel) at a given latitude and longitude on the Earth’s surface. Let’s summarize:
The coordinate system has its center at the solar panel (called the observer). It is traditional to use a left-handed coordinate system where the \(x\)-axis points north, the \(y\)-axis points east, and the \(z\)-axis points toward the zenith point in the direction of the normal vector to the tangent plane. The coordinate system is called the horizontal coordinate system because it is oriented according to the observer’s local horizon.
Any object can be described in the horizontal coordinate system either in Cartesian coordinates \((x,y,z)\) or spherical coordinates \((r,\theta,\phi)\). Here, \( \phi \in [0,2\pi]\) is the azimuth angle measured from the \(x\)-axis towards the \(y\)-axis, \( \theta \in [0,\pi]\) is the zenith angle from the \(z\)-axis, and \(r\) is the radius.
The horizontal coordinate system is fixed to a location on Earth, not the stars or the sun. Objects in the starry sky are imagined to be located on the celestial sphere at a fixed distance, and their position is described by azimuth and zenith angles (radius is thus ignored in this type of description). Over time, the zenith and azimuth angles for an object in the sky (such as the sun) will change as the object appears to move across the sky due to the Earth’s rotation. Since the horizontal coordinate system is defined by the observer’s local horizon, the same object seen from different locations on Earth at the same time will have different values of azimuth and zenith.
We measure the angles either in radians \((\theta, \phi) \in [0, \pi] \times [0, 2\pi]\) or degrees \((\theta,\phi) \in [0, 180^\circ] \times [0, 360^\circ]\). Another commonly used angle is the solar elevation angle (\(\alpha\)), which is complementary to the zenith solar angle (\(\alpha = 90^\circ - \theta\)), measuring the angle from the horizontal plane towards the \(z\)-axis.
Write a Python function
def solar_elevation_angle(theta)that, given \(\theta\) in degrees, calculates \(\alpha\) in degrees.
# WRITE YOUR CODE HERE (if Python is needed)
If the angles are given in radians, we can use solar_elevation_angle together with np.deg2rad and np.rad2deg. Therefore, we don’t need to create functions for both degrees and radians, as we can easily reuse our functions. If, for example, we want to calculate in radians here:
theta_in_rad = np.pi / 3 # zenith angle given in radians
print(np.pi / 2. - theta_in_rad) # elevation angle given in radians
# elevation angle given in radians
np.deg2rad(solar_elevation_angle(np.rad2deg(theta_in_rad)))
It is recommended to continue calculations in radians but to present results/angles in degrees if this is more meaningful/descriptive. Also note that we do not need to use solar_elevation_angle in the following, as both zenith and elevation angles will be available.
In the horizontal coordinate system, any object on the celestial sphere is completely determined by the zenith angle (\(\theta\)) and the azimuth angle \(\phi\). As mentioned, the radial coordinate is ignored when all objects are placed on the celestial sphere. However, we can still include the radial coordinate:
Suppose the sun has a fixed distance \(r_s\) to the Earth. Find a reasonable value for \(r_s\). Provide a (mathematical) expression for how the sun’s \(xyz\)-coordinates can be calculated from \(r_s\), \(\theta_s\), and \(\phi_s\), where \(\theta_s\) and \(\phi_s\) are respectively the zenith and azimuth angles for the sun’s position.
A (flat) solar panel is placed at the origin of the coordinate system. The unit normal \(\pmb{u}_p \in \mathbb{R}^3\) to the solar panel has zenith angle \(\theta_p\) and azimuth angle \(\phi_p\). We consider a normalized/unit solar vector \(\pmb{u}_{s} \in \mathbb{R}^3\) given by \((\theta_s, \phi_s)\). Thus, the sun’s vector field is given by \(\pmb{V} = S_0 \pmb{u}_{s}\).
Provide a (mathematical) expression for \(\pmb{u}_p\) and for \(\langle \pmb{u}_{s}, \pmb{u}_p \rangle\) based on the zenith and azimuth angles. You should simplify the expression so that it contains \(\cos(\phi_p-\phi_s)\) and only 5 trigonometric functions. Show that \(-1 \le \langle \pmb{u}_{s}, \pmb{u}_p \rangle \le 1\). Explain in your own words what it means when \(\langle \pmb{u}_{s}, \pmb{u}_p \rangle < 0\).
Write a Python function
def solar_panel_projection(theta_sun, phi_sun, theta_panel, phi_panel)that returns \(\langle \pmb{u}_{s}, \pmb{u}_p \rangle\) when it is positive and otherwise returns zero.
Take another look at your Python function
def solar_panel_projection(theta_sun, phi_sun, theta_panel, phi_panel). Rewrite it so that it works on NumPy arrays of zenith and azimuth angles. You can test it on the following three situations, where the projection should yield \(0.707107\), \(0.0\), and \(0.0\) (or rather, with numerical errors, it should givearray([7.07106781e-01, 6.12323400e-17, 0.0])). Explain the orientation of the solar panel and the position of the sun in the three situations.
theta_sun = np.array([np.pi / 4, np.pi / 2, 0.0])
phi_sun = np.array([np.pi, np.pi / 2, 0.0])
theta_panel = np.array([0.0, np.pi / 2, np.pi])
phi_panel = np.array([np.pi, 0.0, 0.0])
# solar_panel_projection(theta_sun, phi_sun, theta_panel, phi_panel)
Solar Position Modelling via pvlib#
In Python, the solar position angles, denoted as \((\theta_s, \phi_s)\), can easily be calculated at any location using the Solar Position Algorithm (SPA) with the pvlib package, which is implemented by default with the National Renewable Energy Laboratory’s SPA algorithm [Reda and Andreas, 2003, https://www.nrel.gov/docs/fy08osti/34302.pdf]. We follow https://assessingsolar.org/notebooks/solar_position.html.
import pandas as pd
import pvlib
from pvlib.location import Location
We first need to define the observer’s/panel’s geographical location. This is done using the object pvlib.location.Location in the library pvlib, where we need to specify latitude, longitude, time zone, and altitude, among other parameters. For simulation, data for \((\theta_s, \phi_s)\) from, for example, the calendar year 2024 is used, but here in this initial exercise, we’ll suffice with data for the current month, April 2024.
timezone = "Europe/Copenhagen"
start_date = "2024-04-01"
end_date = "2024-04-30"
delta_time = "Min" # "Min", "H",
# Definition of Location object. Coordinates and elevation of Amager, Copenhagen (Denmark)
site = Location(
55.660439, 12.604980, timezone, 10, "Amager (DK)"
) # latitude, longitude, time_zone, altitude, name
# Definition of a time range of simulation
times = pd.date_range(
start_date + " 00:00:00", end_date + " 23:59:00", inclusive="left", freq=delta_time, tz=timezone
# You might need to use closed='right' instead of inclusive='left' in some versions of pandas
)
Choose a location for your solar panel, for example, DTU. Update the above GPS coordinates (measured in DecimalDegrees), altitude, and name to match the chosen location.
We can now determine the solar position based on the horizontal coordinate system location in site for the specified time interval with the following call:
# Estimate Solar Position with the 'Location' object
sunpos = site.get_solarposition(times)
# Visualize the resulting DataFrame
sunpos.head()
| apparent_zenith | zenith | apparent_elevation | elevation | azimuth | equation_of_time | |
|---|---|---|---|---|---|---|
| 2024-04-01 00:00:00+02:00 | 117.855212 | 117.855212 | -27.855212 | -27.855212 | 339.197342 | -3.866396 |
| 2024-04-01 00:01:00+02:00 | 117.904726 | 117.904726 | -27.904726 | -27.904726 | 339.473686 | -3.866190 |
| 2024-04-01 00:02:00+02:00 | 117.953602 | 117.953602 | -27.953602 | -27.953602 | 339.750301 | -3.865984 |
| 2024-04-01 00:03:00+02:00 | 118.001838 | 118.001838 | -28.001838 | -28.001838 | 340.027184 | -3.865777 |
| 2024-04-01 00:04:00+02:00 | 118.049434 | 118.049434 | -28.049434 | -28.049434 | 340.304332 | -3.865571 |
We see that the DataFrame contains the solar position for each minute in April 2024. The time sampling \(\Delta t\) can be controlled by delta_time = "Min" (minute) set above. When we later need to calculate the energy production for the entire year 2024, it may be sufficient to know the solar position for each hour (for the whole year 2024), as the DataFrame will become very large if we use delta_time = "Min" for a whole year. This is handled by delta_time = "H" (hour). Note that delta_time = "M" (month) sets \(\Delta t\) to a month (which is too large for our needs).
After the solar angles are estimated using pvlib, they can be visualized, for example, for April 1st:
import matplotlib.dates as mdates
chosen_date = "2024-04-01"
# Plots for solar zenith and solar azimuth angles
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(30, 10))
fig.suptitle("Solar Position Estimation in " + site.name + chosen_date)
# plot for solar zenith angle
ax1.plot(sunpos.loc[chosen_date].zenith)
ax1.set_ylabel("Solar zenith angle (degree)")
ax1.set_xlabel("Time (hour)")
ax1.xaxis.set_major_formatter(mdates.DateFormatter("%H"))
# plot for solar azimuth angle
ax2.plot(sunpos.loc[chosen_date].azimuth)
ax2.set_ylabel("Solar azimuth angle (degree)")
ax2.set_xlabel("Time (hour)")
ax2.xaxis.set_major_formatter(mdates.DateFormatter("%H"))
Note that the \(x\)-axis uses the UTC time zone, so you need to add \(+2\) to get Danish time. The plotted vectors can be printed:
chosen_date = "2024-04-01"
print(sunpos.loc[chosen_date].zenith)
print(sunpos.loc[chosen_date].elevation)
print(sunpos.loc[chosen_date].azimuth)
2024-04-01 00:00:00+02:00 117.855212
2024-04-01 00:01:00+02:00 117.904726
2024-04-01 00:02:00+02:00 117.953602
2024-04-01 00:03:00+02:00 118.001838
2024-04-01 00:04:00+02:00 118.049434
...
2024-04-01 23:55:00+02:00 117.237633
2024-04-01 23:56:00+02:00 117.289913
2024-04-01 23:57:00+02:00 117.341566
2024-04-01 23:58:00+02:00 117.392590
2024-04-01 23:59:00+02:00 117.442983
Freq: min, Name: zenith, Length: 1440, dtype: float64
2024-04-01 00:00:00+02:00 -27.855212
2024-04-01 00:01:00+02:00 -27.904726
2024-04-01 00:02:00+02:00 -27.953602
2024-04-01 00:03:00+02:00 -28.001838
2024-04-01 00:04:00+02:00 -28.049434
...
2024-04-01 23:55:00+02:00 -27.237633
2024-04-01 23:56:00+02:00 -27.289913
2024-04-01 23:57:00+02:00 -27.341566
2024-04-01 23:58:00+02:00 -27.392590
2024-04-01 23:59:00+02:00 -27.442983
Freq: min, Name: elevation, Length: 1440, dtype: float64
2024-04-01 00:00:00+02:00 339.197342
2024-04-01 00:01:00+02:00 339.473686
2024-04-01 00:02:00+02:00 339.750301
2024-04-01 00:03:00+02:00 340.027184
2024-04-01 00:04:00+02:00 340.304332
...
2024-04-01 23:55:00+02:00 337.993232
2024-04-01 23:56:00+02:00 338.267217
2024-04-01 23:57:00+02:00 338.541481
2024-04-01 23:58:00+02:00 338.816021
2024-04-01 23:59:00+02:00 339.090835
Freq: min, Name: azimuth, Length: 1440, dtype: float64
Here we have chosen to plot the solar angles for April 1st.
Plot the zenith, azimuth, and elevation angles of the sun, i.e., \(\theta_s, \phi_s, \alpha_s\), for the entire day of April 20, 2024, as a function of time.
Recommendation: In the following, we will need to work with vectors of, for example, zenith and azimuth angles, especially being able to find maximum values, zero crossings, integrate, etc. Therefore, it is recommended to work with the data from sunpos as NumPy arrays. This can be done by, for example:
np.array(sunpos.loc[chosen_date].elevation)
array([-27.85521243, -27.90472593, -27.95360189, ..., -27.34156572,
-27.39258993, -27.44298348])
Plot the elevation angle of the sun and determine when the sun is highest during the day on April 20, 2024. Explain what it means when \(\alpha_s < 0\) or \(\theta_s > 90^\circ\).
Find the time of sunrise and sunset at DTU on April 20, 2024. Compare with “known” values such as those from DMI. Hint: If you want precise values, you should use
apparent_elevation(apparent sun elevation accounting for atmospheric refraction) instead ofelevation. You do not need to account for the curvature of the Earth.
Find the highest point of the sun in the sky (in degrees) on the summer solstice at DTU, and when during the day it occurs? Hint: You will need to change the start and end dates for the
sunposobject.
Create a Python function that can calculate the highest point of the sun \(\alpha_{max}\) in the sky (in degrees) on a given date (year-month-day) at a given location (e.g., city) specified by latitude and longitude. Hint: The answer should not depend on longitude, as the highest point of the sun in the sky depends only on latitude.
In a previous task, you found an expression for the sun’s \(xyz\) coordinates from \(r_s\), \(\theta_s\), and \(\phi_s\).
Write a Python function (for use with NumPy arrays) that converts from the sun’s zenith and azimuth to the sun’s position given in \(xyz\) coordinates. Remember whether you are working in radians or degrees. The
np.deg2rad()function may be useful. It’s fine to use an approximate value for \(r_{s}\), but you can find a more accurate value with:pvlib.solarposition.nrel_earthsun_distance(times) * 149597870700, where149597870700is the number of meters in an astronomical unit (AU).
Write a Python function that converts from the sun’s position in the sky to zenith and azimuth (in degrees or radians) in \(xyz\) coordinates. The
np.arctan2(y, x)andnp.rad2deg()functions may be useful.
Mathematical considerations#
Flux expressed solely in terms of \(A_0\), \(S_0\), and angles#
In our model, it is assumed that the instantaneous solar irradiance received by a solar panel depends on an atmospheric attenuation factor \(A_0\) (here \(A_0 = 0.5\)) and the maximum irradiance \(S_0\) (e.g. \(S_0 = 1100\,W/m^2\)). For each energy calculation, the orientation of the solar panel is fixed, so that the tilt angle \(\theta_p\) and the azimuth angle \(\phi_p\) are treated as constant variables. The position of the sun at any time \(t\) is given by the zenith angle \(\theta_s(t)\) and azimuth angle \(\phi_s(t)\).
Derive an explicit formula for the instantaneous flux \(\text{Flux}(t,\theta_p,\phi_p)\) received by the panel, where the flux is expressed solely in terms of the following variables and constants:
\(A_0\): the atmospheric attenuation factor (0.5),
\(S_0\): the maximum solar irradiance (e.g. 1100 W/m²),
\(\theta_s(t)\) and \(\phi_s(t)\): the sun’s zenith and azimuth angles at time \(t\),
\(\theta_p\) and \(\phi_p\): the panel’s orientation angles.
\(L\) and \(B\): the length and width of the solar panel
This task builds directly on previous exercises and should be straightforward (without calculations). Remember that the flux must be defined to include only positive contributions, i.e., it must be set to zero (e.g., using a function such as \(\max(x,0)\)) when the sun’s rays hit the back side of the panel. Also, the flux should be set to zero if the sun has not risen — for this, “curly bracket”-notation (“definition by cases”) and a condition on \(\theta_s(t)\) can be used.
Continuous solar angle functions of time#
We assume that the sun’s zenith angle \(\theta_s(t)\) is a continuous function of time \(t\).
Discuss why this assumption is reasonable. (A formal proof from mathematical physics is not expected, but a discussion based on intuitive understanding of the physical model.)
The sun’s azimuth angle \(\phi_s(t)\) is likewise a continuous function of time \(t\), except for the jump from 360 degrees to 0 degrees.
Explain why we can therefore conclude that \(\cos(\phi_s(t)-\phi_p)\) and \(\sin(\phi_s(t))\) are continuous functions of time.
Discontinuities in the flux (over a single day)#
Argue that the solar panel’s flux \( ext{Flux}(t,\theta_p,\phi_p)\) as a function of \(t\) is piecewise continuous. Explain why the flux can have discontinuities. Describe a case where the flux has no discontinuities (over a day) and a case where the flux has exactly one discontinuity (over a day). State how many discontinuities the flux can have at most over one day. If we are in a case where the flux has no discontinuities, does it then automatically follow that the flux is also differentiable?
Approximation of solar angles over a single day#
We consider the sun’s path for a single day — here chosen as April 20, 2024 — and assume that the sun’s zenith angle can be approximated by a trigonometric function of the form
and where:
\(M_z\) is the average zenith angle (mean value),
\(A_z\) is the amplitude (i.e., half the variation around \(M_z\)),
\(\omega_z\) is the angular frequency. If \(t\) is given in hours, then for a 24-hour period one can choose \(\omega_z \approx \frac{2\pi}{24}\),
\(T_z\) is the time shift that places the minimum zenith angle (sun’s highest position) correctly, typically around noon.
For the sun’s azimuth angle we assume a linear (affine) function \(\widetilde{\phi}_s(t) := a\,t + b\) as an approximation to \(\phi_s(t)\), where \(a\) and \(b\) are chosen such that the model satisfies:
Use solar position data for April 20, 2024 (retrieved using pvlib) to determine approximate values for the parameters \(M_z\), \(A_z\), \(\omega_z\), \(T_z\), \(a\), and \(b\). To do this, you must fit the respective functions to the solar position data. Finally, create two separate plots that illustrate your approximation:
One for the zenith angle
\[\begin{equation*} \theta_s(t) \approx M_z + A_z\,\cos\Bigl(\omega_z\,(t - T_z)\Bigr), \end{equation*}\]and one for the azimuth angle
\[\begin{equation*} \phi_s(t) \approx a\,t + b. \end{equation*}\]
(You are welcome to propose an improved model \(a\,t + b\) for \(\phi_s(t)\))
Time integration of the flux over a single day#
We define the total energy absorbed by the panel over a day as
where \(t_1\) and \(t_2\) mark the start and end of April 20, 2024, and \(\eta\) is the panel’s efficiency.
We wish to maximize the energy production \(E(\theta_p, \phi_p)\). We assume the panel faces south so that \(\phi_p = \pi = 180^\circ\), and therefore only need to maximize \(E(\theta_p, \pi)\) as a function of \(\theta_p\). Ideally, we want to maximize the energy production over an entire year, so \(t_1\) and \(t_2\) would mark the start and end of the year, but in this exercise we only consider April 20, 2024. Furthermore, we will replace \(\theta_s(t)\) and \(\phi_s(t)\) with the approximations \(\widetilde{\theta}_s(t)\) and \(\widetilde{\phi}_s(t)\) found in the previous task.
In this subtask we assume (for simplicity) that \(\eta\,L\,B =1\). Choose \(t_1\) and \(t_2\) as the times when \(\widetilde{\theta}_s(t)=\pi/2=90^\circ\) (i.e., approximate sunrise and sunset).
Under these assumptions, derive an expression for \(E(\theta_p, \pi)\). Determine whether \(E(\theta_p, \pi)\) is a continuous and/or differentiable function with respect to \(\theta_p\). Find the maximum of \(E(\theta_p, \pi)\) for April 20, 2024 and state the angle \(\theta_p\) that maximizes the energy production.
Note: In this project, we aim to find the optimal solar panel angle \(\theta_p\) that maximizes energy yield. In Mathematics 1b, we typically perform maximization analyses by identifying stationary points and any exceptional cases. Therefore, it is essential to determine whether \(E(\theta_p, \pi)\) is differentiable with respect to \(\theta_p\), as this differentiability enables precise identification of stationary points and potential exceptions.
Power and Energy Calculations#
Calculating exact integrals can be challenging, especially when the integration needs to be performed over an entire year. Therefore, we now transition to numerical computations in Python to estimate the energy output. As a result, we will no longer need the approximated functions \(\widetilde{\theta}_s(t)\) and \(\widetilde{\phi}_s(t)\) from the previous two questions.
Recommendation: It is recommended to work with everything in radians. Remember that you can use np.deg2rad or np.rad2deg. If you have angles in degrees, you should use np.rad2deg.
We consider April 20th. sunpos contains solar position data for every minute throughout the day. Since there are 1440 minutes in a day, solar position angles over this day are described by a vector of this length, for example:
sunpos.loc[chosen_date].zenith
2024-04-01 00:00:00+02:00 117.855212
2024-04-01 00:01:00+02:00 117.904726
2024-04-01 00:02:00+02:00 117.953602
2024-04-01 00:03:00+02:00 118.001838
2024-04-01 00:04:00+02:00 118.049434
...
2024-04-01 23:55:00+02:00 117.237633
2024-04-01 23:56:00+02:00 117.289913
2024-04-01 23:57:00+02:00 117.341566
2024-04-01 23:58:00+02:00 117.392590
2024-04-01 23:59:00+02:00 117.442983
Freq: min, Name: zenith, Length: 1440, dtype: float64
We only consider \(\theta_p \in [0, \pi/2]\), as \(\theta_p > \pi/2\) corresponds to tilting the panel with the back facing upwards.
Create a Python function that can calculate the flux of the sun’s vector field through the solar panel’s surface for each minute throughout the day. You should use
solar_panel_projection(theta_sun, phi_sun, theta_panel, phi_panel). Remember to only include solar zenith angles \(\theta_s \in [0, \pi/2]\) (why?) so that the panel’s flux is zero if the \(\theta_s\) values (in a vector likesunpos.loc[chosen_date].zenith) are above \(\pi/2\), i.e., 90 degrees.
To find the energy production from the solar panel, we need to integrate the flux (i.e., the power) over the considered time period. We always work in SI units, but you should provide final results (such as the total energy production) in relevant units (for example, both in joules and \(\mathrm{kWh}\)). For integration, we can use the Trapezoidal rule known from Mathematics 1b. You can use your own implementation of the Trapezoidal rule if you prefer, but here we choose to use Simpson’s rule (https://en.wikipedia.org/wiki/Simpson’s_rule) from the SciPy package. The energy production can be specified per \(\mathrm m^2\) of the panel. Remember to include the efficiency of the solar panel (regarding the flux) as described in the standard assumptions.
from scipy import integrate
flux = np.array(...) # From the previous task
# Remember to take into account the efficiency of the panel, according to the standard assumptions.
# dx=60 since there are 60 s between time samples
integral_value = integrate.simpson(..., dx=60)
integral_value
# WRITE YOUR CODE HERE (if Python is needed)
flux here is the vector (np.array) containing the flux calculated for each minute throughout the day. The parameter dx=60 tells SciPy that the flux is sampled every minute (1 minute is 60 seconds in SI unit). If you later choose to use delta_time = "H" when calculating the energy production for a whole year, remember to inform integrate.simps about the changed time sampling, namely dx=3600.
Point the solar panel towards the south, i.e., azimuth angle \(\phi_p = 180^\circ\). Calculate the energy production for April 20th for each integer angle \(\theta_p\) between 0 and 90 degrees.
Optimal Angle#
Now, let’s finally consider the energy production for the entire year 2024. Call a new sunpos object with the relevant time interval.
Point the solar panel towards the south, i.e., azimuth angle \(\phi_p = 180^\circ\). Calculate the energy production for the entire year 2024 for each integer angle \(\theta_p\) between 0 and 90 degrees.
Find the optimal angle \(\theta_p\) and indicate the energy production. How much less is the energy production if \(\phi_p\) is, for example, \(175^\circ\) or similar?
Set up a realistic configuration of \(X\) number of solar panels, where you choose \(X\) according to a typical setup on a single-family house. Solar panels are set up at the optimal angle. Calculate the energy production for each day and plot this as a function of time (specified in days).
Extensions#
You must now choose to work further with at least one of the following two extensions:
Solar Panel Film on a Curved Surface#
Solar panel film (thin solar/power film) is a thin film that can be attached to various building surfaces and functions as a solar panel. Go on an excursion in the vicinity of DTU, find a non-planar surface in the urban environment suitable for mounting solar panel film. By a “non-planar” surface, we mean a curved surface where the surface’s normal vector is not constant. It could be a roof of a bus shelter, the top floor of the Jægersborg water tower, etc. It should be a surface that you can parameterize. Set up a parameterization for the chosen building surface. Find a datasheet for a solar panel film that can reasonably be used on the selected building surface and provide the relevant data. Generalize your model and Python code to account for the fact that the solar panel’s normal vector is not constant. Calculate the annual energy production for the solar panel film on the building component.
Since the normal vectors are no longer constant over the entire surface, you must now integrate the projection of \(\pmb{V}\) onto the surface’s unit normal vector over the entire surface. For this, you can use SciPy. Remember that the projection should be set to zero when it is negative. You should choose a surface where the shadow for the surface itself does not become too complicated to handle.
Optimization with respect to profit#
Instead of solely maximizing the annual energy production (in kWh) for the solar panels, you must now optimize the panel orientation (both \(\theta_p\) and \(\phi_p\)) in order to maximize the annual profit, assuming that all electricity is sold to the Danish power grid. Profit is calculated as the total revenue, where the energy produced in each time interval is multiplied by the time-dependent electricity price (in DKK/kWh). Electricity prices vary hourly and can be retrieved via the API from Energidataservice.dk.
The profit function \(P(\theta_p,\phi_p)\) can be formulated as:
where \(\text{Flux}(t,\theta_p,\phi_p)\) (in W) is, for example, converted into kWh per time interval if \(\text{Price}(t)\) is the electricity price in DKK/kWh for the given time interval.
Derive and/or implement the profit function \(P(\theta_p,\phi_p)\) for an entire calendar year (e.g., using data from 2024), where electricity prices vary hourly.
Determine the optimal panel angle \(\theta_p\) and \(\phi_p\) that maximize the profit \(P\).
Discuss how the time-varying price structure may change the optimal orientation compared to optimization based solely on energy production. For instance, whether it might be advantageous to adjust the panel’s azimuth angle to produce more energy during hours with high electricity prices.
Assumptions and considerations:
Use the same models for flux and solar position calculations as before.
Electricity price data can be retrieved via the API at Energidataservice.dk.
We assume that the electricity can be sold without additional costs (fees, transport, etc. are ignored) for a price \(\text{Price}(t)\) that depends linearly on the spot price (for DK1).
Short guide for retrieving price data:
Check the documentation on Energidataservice.dk to find the endpoint for hourly prices.
In Python, you can use the
requestslibrary to send a GET request. A simple example might be:
import requests
import pandas as pd
endpoint = "https://api.energidataservice.dk/dataset/elspotprices"
params = {
"start": "2024-01-01T00:00",
"end": "2024-12-31T23:00",
"filter": '{"PriceArea":["DK2"]}' # <-- fixed JSON string
}
response = requests.get(endpoint, params=params)
data = response.json()
df = pd.DataFrame(data['records'])
# Convert units if needed (e.g., from øre to DKK, or J/MJ to kWh)
print(df.head(20))
Integrate the retrieved hourly prices into the profit function \(\text{Price}(t)\).
Optional extensions#
Include atmospheric attenuation#
The current model assumes that the solar irradiance \(S_0\) is constant (e.g., 1100 W/m²) for all times when the sun is above the horizon. This means that if a solar panel is oriented such that \(\pmb{u}_p\) is parallel with the sun’s rays, it will receive the same flux—regardless of the time of day, e.g., sunrise or noon. But in reality, this is not the case:
Atmospheric attenuation: Especially at sunrise and sunset, the sun’s rays must pass through a longer stretch of the atmosphere, leading to increased scattering and absorption. This reduces the actual irradiance received by the panel.
Temperature and spectral distribution: The sun’s radiation may change character depending on the time of day—for instance, the shortwave part of the spectrum is more intense around midday, which affects panel efficiency.
Air mass effect: The so-called “air mass” is a measure of how much atmosphere the sun’s rays must pass through. At low solar elevations, the air mass is greater, and therefore the effective irradiance is lower than it is at high solar elevations (e.g., around noon).
To make the model more realistic, consider including an attenuation factor that depends on solar elevation (or air mass). This means that orientation toward the sun at sunrise will not yield the same flux as at noon, even if the panel orientation is “optimal” in the sense of being directly hit by the sun’s rays.
Solar energy during colonization of Mars and Europa (Jupiter’s moon)#
Assume we want to optimize the output of solar panels during a potential colonization of Mars—and then examine the feasibility of solar energy during colonization of Europa (Jupiter’s moon). In the current model, the solar irradiance \(S_0\) is treated as a constant (e.g., 1100 W/m²) for Earth. How should \(S_0\) be adjusted for Mars and Europa based on the inverse-square law, given that, e.g., Mars’ average distance to the sun is approximately 1.524 times Earth’s distance? Answer the following:
For Mars:
Compute the new solar irradiance \(S_{\text{Mars}}\) on Mars using the inverse-square law, given that Mars’ average distance to the sun is about 1.524 times Earth’s distance:
\[\begin{equation*} S_{\text{Mars}} = S_{\text{Earth}} \times \left(\frac{d_{\text{Earth}}}{d_{\text{Mars}}}\right)^2. \end{equation*}\]Discuss how the lower irradiance affects total energy production from solar panels on Mars, and what additional factors (e.g., Mars’ thin atmosphere, dust storms, and temperature variations) must be considered.
For Europa (Jupiter’s Moon):
Calculate the expected solar irradiance \(S_{\text{Europa}}\) for Europa, noting that the average distance to the sun for Jupiter’s system is about 5.2 times Earth’s distance:
\[\begin{equation*} S_{\text{Europa}} = S_{\text{Earth}} \times \left(\frac{d_{\text{Earth}}}{d_{\text{Europa}}}\right)^2. \end{equation*}\]Discuss whether the calculated irradiance on Europa provides a realistic foundation for solar energy utilization. Consider challenges such as lack of a dense atmosphere, surface characteristics, and high radiation exposure from Jupiter.
Summary:
Compare the computed values of \(S_{\text{Mars}}\) and \(S_{\text{Europa}}\) to Earth’s irradiance \(S_{\text{Earth}}\) (e.g., 1100 W/m² under ideal conditions).
Discuss how the reduction in solar irradiance affects the design and placement of solar panels in the two colonization scenarios, and whether solar energy can be a viable energy source under these conditions.
Hint: Use the following formula for both calculations:
where \(d_{\text{new}}\) is the average distance to the sun for the respective celestial body (1 AU for Earth, 1.524 AU for Mars, and approx. 5.2 AU for Europa).
Include Shadow from a Tree or Building#
In the horizontal coordinate system, place, for example, a tree crown at a given distance, such as \(10 \,\mathrm m\). Consider the sphere \(r_{tree}=10\,\mathrm m\) with the center at the Origin (the panel’s position) and describe the tree’s shape on this sphere. A simple model is as follows: Assume that there are only two possibilities: either \(0\%\) shadow or \(100 \%\) shadow. Calculate the surface area of the entire sphere as a function of \(r_{tree}\), assume a size of the tree crown (e.g., its approximate diameter), and calculate approximately where the tree shades in zenith and azimuth degrees. This could be that there is \(100 \%\) shadow when \(\theta \in [70,80]\) and \(\phi \in [150,160]\), and the energy intake must therefore be subtracted when the sun is in this interval. The shape of the tree becomes somewhat unnatural when using axis-parallel areas in the \((\theta,\phi)\) plane. What shape does it correspond to on the sphere \(r_{tree}=10\,\mathrm m\)? Discuss how much effect shadows from trees and buildings can have. Is it possible to describe a perfectly spherical tree crown? How would this shape look on the sphere \(r_{tree}\)?
Alternatively, you can model shadow from a neighboring building. The idea is the same, but the shadow interval for \(\theta\) should go all the way down to the ground, i.e., \(\theta = 90\). Depending on the building size and distance to the solar panel, the angle intervals may well be larger, e.g., \(\theta \in [65,90]\). You can find inspiration with this tool: https://www.findmyshadow.com/
More on Optimization#
Is South-Facing Panel Optimal? It should be, but let’s investigate it mathematically. In the above task, we assumed \(\phi_p = 180^\circ\). You should now drop this assumption and calculate the panel’s energy production as a function of both \(\phi_p\) and \(\theta_p\).
A Panel with a Motor: Assume that the solar panel’s angle can be adjusted either daily, monthly, or quarterly. If you choose, for example, monthly adjustments, then the optimal angle for each month should be specified, and the panel’s annual power should then be calculated. Compare the annual power with the power for a corresponding fixed-mounted panel. Discuss whether it is worth having solar panel systems with panels whose angle can be adjusted in this way.