Optimizing solar panel placement

2023-07-02mathoptimizationpythonsolar

I have a solar panel, and I want to figure out where to install it to capture the most energy. Instead of moving the panel to a new site, adjusting its orientation, and waiting a year to collect data, we can use computational models to provide some guidance. In this post, I’ll present a series of models of increasing sophistication that could be used to optimize solar panel orientation and placement.

First orbit: orientation §

A reasonable place to start is to point the solar panel toward the sun—I mean, let’s face it, the panel isn’t going to do what we want if we put it face down on the ground. So pointing at the sun seems good.

But for a computational model, we need a way to map “pointing at the sun” to some numbers, so we can make things quantitative. To do this we can describe the position of the sun in the sky by adopting a coordinate system.

Coordinate system §

Probably the most widely used coordinate system is the Cartesian coordinate system, which uses two distance values to indicate the horizontal and vertical position of a point on a two-dimensional plane with respect to some predefined “origin.” This isn’t ideal for our purposes, because the sky is better modeled as a spherical surface at some fixed distance from an observer. Circles and spheres suggest that angles, instead of distances, will be more useful for comparing multiple locations in the sky.

A diagram of an azimuth-altitude coordinate system.^[1]

A horizontal coordinate system seems to be just the thing then. Because the surface of a sphere is two-dimensional, this coordinate system uses two angular values to describe the location of an object in the sky:

• azimuth – the angle of an object along the horizon. We’ll measure it in degrees, with north being 0°, increasing as one turns toward the east.
• altitude – the angle of an object above the horizon. We’ll measure it in degrees, with the horizon being 0°, increasing as the object rises higher in the sky. (Negative altitude values are possible and indicate that an object is hidden below the horizon.)

As an example, let’s follow the sun over the course of a day. At sunrise, the sun is on the horizon in the east, or in our coordinates, at an altitude of 0° and an azimuth of 90°. Over the course of the morning, the sun travels higher up into the sky, covering altitude values between 0° and some peak, let’s say 50°. At the same time, the sun traverses the southern part of the sky (assuming the observer is in the northern hemisphere), so its azimuth values change between approximately 90° (east) and 180° (south). In the afternoon, the sun’s altitude starts decreasing from 50° back toward the horizon at 0°, as the azimuth values continue to increase from 180° (south) toward sunset around 270° (west).

Of course, the exact values here depend on the observer’s location on the earth, as well as the season of the year. But in general terms, this coordinate system gives us a way of describing the location of the sun in the sky, for any time of day, any time of year, and even any location on earth.

One of the great things about using this coordinate system is that it maps pretty well onto the orientation of our solar panel. For example, let’s suppose it’s 10 in the morning, and the sun is at an azimuth of 120° and an altitude of 20°. We could rotate our solar panel to point at the sun by facing it toward 120°, using something like a compass. Then we need to stand the panel up (pointing it at the horizon), and lean it back against a wall or something, so that it’s tilted 20° up from the horizon.

Movement §

An example of a page from an ephemeris.^[2]

This all sounds exciting so far, but there’s a problem hinted at by the sunrise-to-sunset example: the sun moves in the sky. This is a problem for two reasons actually:

1. We need to know how the sun moves through the sky over the course of a day.
2. Because the sun is moving, we need to decide how to orient the solar panel so that it will capture as much power as possible over time.

It turns out that the first problem (how the sun moves) is solved already. The location of the sun in the sky has been important to humans since… forever, so the paths of celestial bodies through the sky have been documented for millenia; a book or table of these quantities is known as an ephemeris. There are many historical and contemporary ephemera, like books and software programs, describing the azimuth and altitude of the sun at any given moment. So, thanks to our fellow humans, we can compute the location of the sun in the sky for a whole year, without having to put the solar panel outside and wait for the sun to travel around.

The second problem (that the sun moves) has also been solved, actually, in a rather obvious way: just point the solar panel directly at the sun at every moment! But doing this requires us to go out and move the solar panel around all the time, which is generally annoying. An alternative is to build a machine that will move the solar panel around for us. Such a machine is called a solar tracker, which does exist already, but tends to be expensive and introduce complexity—as well as the possibility of mechanical failures—into our solar setup. So, for the moment, I’m going to say that I want to work out a model to identify a fixed orientation for my panel. But I’ll get back to trackers later on, once we have a good model in place.

Model §

So if my panel is going to be in a fixed orientation, but the sun is moving, then we need a way to measure how well the panel is capturing power from the sun, when the panel is only pointed “near” the sun. Time for some math!

Let’s start out by considing an “overhead” view of the panel and the sun, so we can focus just on the azimuth. Our panel will point in a direction, say something like southeast or west, and we’d like to measure how closely the panel is pointed toward the sun when the sun is, say, at an azimuth of south or southwest.

Intuitively, if the sun is at 180° (south), then we’d probably want to point the panel also at 180°. But what if we couldn’t get it to be exactly 180°? How much would we lose by pointing the panel toward 170° or 30°?

Hopefully it makes sense that if we pointed the panel northward (0°), it wouldn’t get any direct sunlight when the sun was in the south (180°). Hopefully it also makes sense that the panel still wouldn’t get any direct light if we pointed it toward the east (90°) or west (270°), since the face of the panel would be parallel to the sun’s rays at that point. But anywhere between east and south would be able to catch some light, probably with better performance the closer we got toward south.

A common way to model this relationship is to use the cosine function of the angle between two things, as a way to quantify how closely they are pointed in the same direction. Let’s represent the solar azimuth at time \( t \) as \( \alpha_☼(t) \) and the panel azimuth as \( \alpha_▱ \). Using a cosine distance, the alignment \( \xi \) of the panel with the sun is given by \[ \xi(t) = \cos( \alpha_☼(t) - \alpha_▱ ). \] The cosine function is great for this for a number of reasons: it is defined over angles, it varies smoothly, it’s symmetric around zero, and its output ranges from zero (when the panel is perpindicular to the sun) to one (when the panel is pointed directly at the sun).

Two dimensions §

Recall that our horizontal coordinate system uses not one but two numbers to represent the position of the sun in the sky: the azimuth \( \alpha_☼(t) \) as well as the altitude, which we’ll represent as \( \theta_☼(t) \). I’ll also introduce \( \theta_▱ \) for the altitude of the panel, so that the panel altitudes match those of the sun, i.e., a panel altitude of 0° means that the panel is pointing at the horizon, while an altitude of 90° means the panel is flat on the ground, pointing straight up at the sky.

One of the nice things about the cosine function for measuring angular similarity is that it decomposes nicely in two dimensions. So if we measure the azimuth similarity as above, using the cosine of the angular difference in azimuth values, then we can also measure the altitude similarity using cosine: \[ \cos( \theta_☼(t) - \theta_▱ ). \] To combine the two cosine values, note that

1. we want both of the angles to match as closely as possible, and
2. if either angle is mismatched, then the overall performance will be reduced.

To check our intuition on point 2, imagine that our panel matches the solar azimuth (that is, \( \alpha_☼(t) = \alpha_▱ \) ), but the altitude is 90° off (so, for example, our panel is pointed at the horizon, while the sun is directly overhead). In this case, the panel would receive no light from the sun, even though the azimuth matches perfectly!

So we want a combination that maps to zero if either individual cosine value is zero, and maps to one only when both individual cosine values are one. This combination can then be achieved by multiplying together the two cosine metrics: \[ \xi(t) = \cos(\alpha_☼(t) - \alpha_▱) \cos(\theta_☼(t) - \theta_▱). \]

Integration §

Instead of maximizing the performance of my panel for one instant, I’d like to optimize the performance over a whole time period—like, say, a specific month, or an entire year. Again assuming the azimuth and altitude of the panel are fixed, while the azimuth and altitude of the sun vary with time, we can define an integration of the panel output from starting time \( t_i \) to ending time \( t_f \) as \[ E_▱ = \sum_{t_i,t_i+\Delta t,\dots,t_f} \xi(t) \] where \( \Delta t \) is some small-ish time increment.

The task that I’d like the computer to handle is to figure out which values of \( \alpha_▱ \) and \( \theta_▱ \) give the maximum output for \( E_▱ \).

Coding it up §

I think we have the ingredients that we need to enlist a computer’s help with our problem!

Remember that I mentioned an ephemeris? We need one to compute the position of the sun in the sky. There are several Python packages that do this, but I’ll use one called suncalc. We need to provide suncalc.get_position with the time that we’d like to query, as well as the panel’s location on the surface of the earth.

import collections
import numpy as np
import suncalc

LATITUDE, LONGITUDE = 37, -122

Coords = collections.namedtuple(
    'Coords', 'azimuth_deg altitude_deg')

def sun_position(t):
    r2d = np.rad2deg
    sun = suncalc.get_position(t, LONGITUDE, LATITUDE)
    azimuth = (np.pi - sun['azimuth']) % (2 * np.pi)
    return Coords(r2d(azimuth), r2d(sun['altitude']))

Note that suncalc.get_position takes the longitude value first, and the latitude second! I’m also converting the azimuth value returned by suncalc so that it matches the north-to-east convention that I described at the start of this post.

Next, we can define a function that computes the alignment \( \xi \) between the sun and the solar panel (or, in general, any two points in the horizontal coordinate system):

def alignment(u, v):
    c = lambda a, b: np.cos(np.deg2rad(a - b))
    azimuth = c(u.azimuth_deg, v.azimuth_deg)
    altitude = c(u.altitude_deg, v.altitude_deg)
    return azimuth * altitude

We can integrate from a start to an end time, as a starting point, by summing up the alignments between the sun and the panel at each moment:

import datetime

DT_H = 0.5  # Time step in hours.

def summed_alignment(panel, t_start, t_end):
    acc = 0
    t = t_start
    while t < t_end:
        sun = sun_position(t)
        t += datetime.timedelta(hours=DT_H)
        if sun.altitude_deg < 0:
            continue
        acc += alignment(panel, sun)
    return acc

Note that the time that we provide to suncalc.get_position needs to be in UTC, not local time, so we’ll need to convert a local time value into UTC using the zoneinfo library.

import zoneinfo

TIMEZONE = zoneinfo.ZoneInfo('America/Los_Angeles')
UTC = datetime.timezone.utc

def utcdate(y, m=1, d=1):
    return datetime.datetime(
        y, m, d, 0, 0, 0, tzinfo=TIMEZONE
    ).astimezone(UTC)

Finally, we can invoke the minimizer:

import scipy.optimize

if __name__ == '__main__':
    opt = scipy.optimize.minimize(
        lambda x, *args: -summed_alignment(Coords(*x), *args),
        args=(utcdate(2019), utcdate(2020)),
        x0=(120, 20),
        bounds=((0, 360), (0, 90)),
        callback=print,
    )

    print('Optimal solar panel orientation:', opt.x)
    print('Total alignment:', -opt.fun)

We provide the minimize function with the function to minimize, as well as an initial guess for our panel orientation. I’ve also provided bounds to the optimizer, one (lower, upper) tuple per variable that we’re optimizing. This can be handy if there are constraints on the panel’s orientation for whatever reason. Note that we actually minimize the “cost,” here formulated as the negative converted energy.

And, with that, congratulations, we’re done! 🎉 The code above has the basic ingredients for computing an optimal solar panel orientation. Time to harvest some electricity from the sun! 💡

Full script by the end of orbit 1: panel-optimization-orbit1.py

Second orbit: location §

Unless my solar panel is at sea, or I live in a place with no hills and no trees, the model developed in the first orbit has a pretty glaring inaccuracy—I’ve assumed so far that the horizon is at an altitude of 0° at every azimuth value! For the rest of us, there could be some obstacles blocking out some of our sunlight, like:

• Nearby topography. For example, if you live in the bottom of a valley, the mountains on either side could block out the sun at some times of the day.
• Nearby tall objects. Trees, tall buildings, or other above-ground features could block out the sun at some times of the day.

To satisfy the constraints imposed by the real horizon in a particular location, it would be best to find an available place (e.g., on the roof, or in the backyard, or whatever we have to choose from) that doesn’t block out views of the sky.

But not all parts of the sky are equally important to keep clear! For example, in the northern hemisphere, the sun will be in the northern half of the sky much less often than in the southern half—and vice-versa for the southern hemisphere. I live in the northern hemisphere, so if I have a space that backs up to a bunch of tall trees on the north, that’s probably ok, but if, say, my roof is shaded by trees to the south, that could be more difficult.

So what we need to know is, for each point on the horizon, what altitude does the sun need to achieve to shine onto the solar panel? Once we have this information we can modify the model to sum up only the solar energy that the panel can actually receives, when the sun’s altitude at a given azimuth exceeds the (apparent) horizon’s altitude at that azimuth.

Measuring the apparent horizon §

An inclinometer is used to measure angles with respect to the direction of gravity.^[3]

There are a couple of steps to this process. The first is that we need to go outside with some tools and measure the altitude values of the apparent horizon, as seen from the place where the solar panel will go. This only takes a little time, and the following tools:

• pen & paper
• compass, or (somewhat better!) a compass app
• inclinometer, or a level or accelerometer app that can show degrees of tilt

The basic process goes as follows:

1. Use the compass to orient to a particular azimuth value. (If using a magnetic compass, be sure to add the magnetic declination, since we really want to be using the geographic azimuth values here.)
2. Sight along the inclinometer or tilt meter until it is higher than the apparent horizon in this direction.
3. Record the altitude of the apparent horizon at this azimuth.
4. Repeat every N° of azimuth, until you’ve covered all the necessary azimuth values.

An inclinometer can be used to measure the angles of distant objects.^[4]

Fitting a curve §

Now, armed with a set of azimuth and altitude values, we need to bring them into our code. Again, what we want is the apparent altitude of the horizon at a given azimuth. Since we have some finite collection of measurements to span the horizon, we could take the measured data and mathematically connect the dots, to form a piecewise linear approximation of the horizon. But as long as we’re connecting the dots, we might as well use a spline, which will connect the dots smoothly!

import scipy.interpolate

# Azimuth values where I took horizon
# measurements. East to west, every 10°.
AZIMUTHS_DEG = np.arange(90, 271, 10)

def _horizon_spline(*altitudes_deg):
    return scipy.interpolate.UnivariateSpline(
        AZIMUTHS_DEG, altitudes_deg)

SITES = dict(
    a=_horizon_spline(43, 36, 30, 30, 17, 30, 30, 28, 24, 23, 18, 22, 18, 17, 14, 18, 19, 20, 22),
    b=_horizon_spline(10, 10, 10, 8, 9, 8, 10, 10, 16, 17, 38, 38, 29, 29, 23, 25, 26, 43, 43),

    # This "site" is a best-case
    # scenario with a flat horizon.
    flat=_horizon_spline(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),
)

Adding to the code §

In general, we can add the horizon information as an extra argument to most of the functions we defined in orbit 1. The place where the horizon needs to be used is in our integral, on the line

if sun.altitude_deg < 0:

This is where I was making an assumption that the horizon is always at zero! Instead, we need to compare to the true horizon altitude, for the computed sun’s azimuth. To get the horizon to this point in the code, we also need to pass it as an extra argument to the integration function:

def summed_alignment(panel, t_start, t_end, horizon):
    acc = 0
    t = t_start
    while t < t_end:
        sun = sun_position(t)
        t += datetime.timedelta(hours=DT_H)
        if sun.altitude_deg < horizon(sun.azimuth_deg):
            continue
        acc += alignment(panel, sun)
    return acc

And we can loop over the potential available horizons in our main driver:

if __name__ == '__main__':
    for name, horizon in SITES.items():
        opt = scipy.optimize.minimize(
            lambda x, *args: -integral(Coords(*x), *args),
            args=(utcdate(2019), utcdate(2020), horizon),
            x0=(120, 20),
            bounds=((0, 360), (0, 90)),
            callback=print)

        print('Site:', name)
        print('Optimal solar panel orientation:', opt.x)
        print('Total alignment:', -opt.fun)

Full script by the end of orbit 2: panel-optimization-orbit2.py

Third orbit: physical grounding §

I like to ground my models as much as possible in the real world, so it’s time to start considering units: how will I measure the output of my solar panel? I think it’s particularly useful to think about units for this task, because our optimizer is able to output a value, but it’s not really possible to gauge intuitively whether that value is in the right general range or not. Luckily, we’re dealing with a physical system, so we can map our model onto appropriate physical units.

Units §

The sun’s light provides energy to the earth per unit of time and per unit of surface area. Intuitively, these units seem ok: we get more energy from the sun the longer we have a panel out in the sunlight, and we get more energy from the sun if we cover more of the roof or backyard with solar panels.

Energy is measured in Joules (J), time is measured in seconds (s), and we’ll measure surface area in square meters (m²). Power, measured in Watts (W), is defined as energy per unit time, so the energy per unit time on the surface area of a solar panel is often expressed in kilowatts per square meter (kW/m²) instead of 1000 J/s/m². A ballpark figure for this solar power per unit area (called “surface power density”) is \( I \) = 1 kW/m², but the main problem with this estimate is that it is constant—remember that the sun moves through the sky over time, and the actual surface power density available at a given moment depends on the sun’s altitude in the sky.

Atmosphere §

The sun continually radiates energy into space, and some of this energy reaches the earth’s atmosphere after about seven minutes. It incurs some loss during this trip, but we’ll consider the solar power density arriving at earth to be a constant. The energy that reaches the outside of the earth’s atmosphere must then travel through some part of the atmosphere before it reaches the surface of the planet. Traveling through the atmosphere causes additional loss, due to things like scattering or absorption, and the losses will be greater the more atmosphere the light has to travel through. For example, consider how the sun feels less strong at sunrise or sunset than at noon.

The solar energy concept that’s used to model atmospheric energy loss is called air mass, abbreviated AM. Air mass represents the relative increase in atmosphere that sunlight must traverse, relative to sunlight entering from directly overhead (the zenith). According to Wikipedia, there are several models for mapping zenith angle to air mass, but I picked one that seems to balance simplicity and accuracy: \[ AM(z) = \frac{2r+1}{\sqrt{(r \cos z)^2 + 2r+1} + r \cos z} \] where \( r = R_E/y_a \) is the ratio of the radius of the earth (approximately 6371 km) to the thickness of the atmosphere (assumed to be 9 km for the purposes of this model). The zenith angle is the complement of the altitude, i.e., \( z = \frac{\pi}{2} - \theta \).

I’m afraid that describing the model in such detail just makes it seem more complex than it is. Fear not! All of the above translates in a straightforward way to Python code, something like:

EARTH_RADIUS_KM = 6371
ATMOSPHERE_KM = 9

def air_mass(altitude_deg):
    z = np.pi / 2 - np.deg2rad(altitude_deg)
    r = EARTH_RADIUS_KM / ATMOSPHERE_KM
    r_cos_z = r * np.cos(z)
    return (2*r + 1) / (np.sqrt(r_cos_z**2 + 2*r + 1) + r_cos_z)

When the sun is at an altitude of 90°, or a zenith angle of 0°, the air mass is equal to 1; any deviation from the zenith results in an air mass greater than 1. Note: For consistency with the remainder of this post, I’ve included a conversion here from altitude to zenith angle, so that we can compute the air mass using an altitude value as an input.

The next tool we need is to map the air mass, \( AM(\theta_☼(t)) \), onto an incident power level that reaches the surface of the earth, \( I(\theta_☼(t)) \). Wikipedia again comes to the rescue with this model:

\[ I(t) = I(\theta_☼(t)) = I_0 (1 + d) 0.7^{AM(\theta_☼(t))^{0.678}} \]

where \( I_0 \) is the power density from the sun in space (outside the atmosphere) and \( d \) is the relative fraction of diffuse (as opposed to direct) light reaching the surface. Some Python code to implement this could look like:

DIFFUSE_FRACTION = 0.1
SPACE_POWER_DENSITY_KW_M2 = 1.353

def surface_power_density_kw_m2(altitude_deg):
    if altitude_deg < 0:
        return 0
    return (SPACE_POWER_DENSITY_KW_M2 *
            (1 + DIFFUSE_FRACTION) *
            0.7**(air_mass(altitude_deg)**0.678))

The value of the power density in space is, unlike the surface measurement, independent of the position of the sun relative to the earth’s surface, and fairly constant at 1361 W/m². This model provides us with additional accuracy, in particular for times when the sun is near the horizon—at such times, the incident solar density on the surface is quite low, and so we’ll be able to take this into account when computing an optimal panel orientation.

Efficiency §

To summarize the model so far, the surface power density that the solar panel captures can be computed by multiplying, for a given moment in time, the available incident power density on the panel, \( I(t) \), by the alignment of the panel and the sun, \( \xi(t) \).

This power density is captured for each square meter of panel that we have available, so we can multiply by the surface area of the panel, \( A \), to get the total captured power.

One final practical issue is that solar panels can only convert a fraction of the solar energy that they capture into usable electricity. This fraction, called the efficiency, tends to come in around 20%, but the exact efficiency value \( \eta \) depends on the panel.

Integration §

We’ll incorporate these constants into our model, giving the following expression for the converted solar power, as a function of time: \[ P(t) = A \eta I(t) \xi(t). \] To measure the total converted solar energy, we’ll need to integrate this model over time, like we did in orbit 1. However, the units here suggest that our previous integration was missing an ingredient: if we are summing power values, we’ll get power as an output; instead, we need to multiply power by time to get energy. This gives us the following integration: \[ E(t_i, t_f) = A \eta \Delta t \sum_{t_i,t_i+\Delta t,\dots,t_f} I(t) \xi(t). \]

Coding it up §

The last step here is to modify our code to incorporate our (a) expression of surface power density, and (b) units in all of these computations. The only term inside the sum that needs to change is the surface power density, since that’s a function of the sun’s position over time. The remaining constants we can just multiply after doing the sum. All of these changes can happen inside the integration function!

def energy_kwh(panel, t_start, t_end, horizon):
    acc = 0
    t = t_start
    while t < t_end:
        sun = sun_position(t)
        t += datetime.timedelta(hours=DT_H)
        if sun.altitude_deg < horizon(sun.azimuth_deg):
            continue
        kw_m2 = surface_power_density_kw_m2(sun.altitude_deg)
        acc += kw_m2 * alignment(sun, panel)
    return acc * DT_H * PANEL_EFFICIENCY * PANEL_AREA_M2

Full script by the end of orbit 3: panel-optimization-orbit3.py

Fourth orbit: weather §

The orbits of our model have progressively refined the assumptions of the original model, improving accuracy and making the model output more interpretable. Somewhat disappointingly, each orbit has progressively reduced the estimated output from the solar panel, which makes sense in a way because the original model was so optimistic about possible limitations, like assuming a flat horizon.

In this orbit I’ll address another assumption of the model by looking at clouds. So far, I’ve assumed that the sun is shining on the solar panel, as long as the sun’s altitude is higher than the apparent horizon. But in reality, the sky is not always perfectly clear. Like the previous orbits, refining this assumpion in the model is going to reduce the estimated output from the panel, but that’s actually a good thing, because it’s better to have a more accurate estimate of the panel’s performance!

Like the path of the sun through the sky, humans have always been rather obsessed with the weather. Unlike the sun, however, the exact conditions of the sky at a given moment are unpredictable, so the best we can do is look at previous conditions of the sky, and hope that the future looks somewhat like the past.

I used weatherspark.com to look up cloud coverage for the 21st day of each month in my location. This site provides frequencies for each cloud coverage level from 10% (almost completely clear) through 30%, 50%, 70%, and 90% (almost completely cloudy). The observed frequencies of all the cloudiness levels add up to 100; that is, for a given day \( d \), we get a probability distribution \( p_d(☁ = c) \) that specifies how likely it is to be \( c \) percent cloudy on that day.

Sampling §

Because the cloudiness on a given day isn’t certain, but we do have a distribution of how likely it is to be cloudy vs. clear, we can use our data in a couple of different ways.

One way is to draw a sample from the cloudiness distribution each time we want to know how cloudy it’s likely to be. So, for example, if we’re trying to compute the solar power that’s available on April 3rd, we can look up the distribution for that day, \( p_{Apr3}(☁) \), pick a random number between 0 and 100, and compare it to the cloudiness distribution.

For example, let’s say the cloudiness distribution for April 3rd is \( p_{Apr3}(☁ = c) \) = [45, 15, 10, 10, 20], meaning that the probability of 10% cloudiness is 45%, the probability of 30% cloudiness is 15%, and so on up to a 20% chance of 90% cloud coverage. First we’ll compute the cumulative distribution of these values, giving \( P_{Apr3}(☁ \le c) \) = [45, 60, 70, 80, 100]. These values indicate the percent of time that the sky has at most the corresponding level of cloud coverage, i.e., 45% of the time it’s at most 10% cloudy, 60% of the time it’s at most 30% cloudy, and so on.

Then if we pick a random number and get, say, 40, we’d consider our sample to correspond to the 10% cloudiness level. If we get a random number like 75, that would correspond to a cloudiness level of 70%.

This process is a perfectly reasonable way to estimate cloudiness using historical weather data. The downside of this approach is that it adds so-called noise to our solar panel calculations. Each time we estimate the output from our solar panel, we might get a different value because we’ve gotten different guesses for how cloudy it might be.

Expectation §

There’s another way to use the information in our cloudiness distribution: compute a weighted sum of the cloudiness levels, using the distribution probabilities as weights. This is called computing an expected value; it gives us the best estimate of the value that we’ll get when we sample randomly from a distribution. Also, the expected value of a distribution is constant, so we’ll eliminate any noise that a sampling procedure might have added.

Mathematically, an expected value looks something like: \[ \mathbb{E}[c] = \sum_c c p_d(☁=c) \] In code it’s quite similar, but since we have percentages as weights, we’ll divide by the sum of the weights to convert these into probabilities.

CLOUD_COVERAGES = np.array((0.1, 0.3, 0.5, 0.7, 0.9))

def _expected_cloud_coverage(*weights):
    assert sum(weights) == 100
    return (weights * CLOUD_COVERAGES).sum() / sum(weights)

I mentioned earlier that I fetched cloud data for the 21st day of each month. I’ll connect those estimates smoothly using a spline that will map the day of the year onto an expected cloudiness level:

def dayofyear(t):
    return t.toordinal() - utcdate(t.year).toordinal()


EXPECTED_CLOUD_COVERAGE = scipy.interpolate.UnivariateSpline(
    [dayofyear(utcdate(2019, m, 21)) for m in range(1, 13)],
    [_expected_cloud_coverage(30, 7, 10, 10, 43),
     _expected_cloud_coverage(28, 9, 10, 11, 42),
     _expected_cloud_coverage(33, 9, 9, 12, 37),
     _expected_cloud_coverage(39, 10, 9, 10, 32),
     _expected_cloud_coverage(46, 12, 10, 9, 23),
     _expected_cloud_coverage(68, 9, 6, 5, 12),
     _expected_cloud_coverage(79, 8, 5, 3, 5),
     _expected_cloud_coverage(78, 7, 5, 3, 7),
     _expected_cloud_coverage(71, 9, 6, 5, 9),
     _expected_cloud_coverage(53, 9, 7, 9, 22),
     _expected_cloud_coverage(35, 9, 9, 10, 37),
     _expected_cloud_coverage(30, 9, 9, 9, 43)]
)

Adding to the code §

This additional limitation of our panel’s output can, like the other time-varying factors, go inside the body of our integral.

def energy_kwh(panel, t_start, t_end, horizon):
    acc = 0
    t = t_start
    while t < t_end:
        sun = sun_position(t)
        t += datetime.timedelta(hours=DT_H)
        if sun.altitude_deg < horizon(sun.azimuth_deg):
            continue
        kw_m2 = surface_power_density_kw_m2(sun.altitude_deg)
        clear = 1 - EXPECTED_CLOUD_COVERAGE(dayofyear(t))
        acc += kw_m2 * clear * alignment(sun, panel)
    return acc * DT_H * PANEL_EFFICIENCY * PANEL_AREA_M2

Full script by the end of orbit 4: panel-optimization-orbit4.py

Fifth orbit: tracking §

I promised earlier that I’d get back to the question of tracking, so here we (finally) are. The model, after four orbits, reflects the path of the sun in the sky, the appearance of the horizon, effects of the atmosphere, and different weather conditions. So far I’ve been obstinately refusing to add a tracker to my panel, but at what cost? We can modify the model slightly to get a quantative answer to this question.

Remember where we compute the alignment of the solar panel with the position of the sun? A tracker ideally follows the path of the sun in the sky, so that the alignment between the two is always perfect—which our model would provide with an alignment score of 1.

We can add this feature to the code by forcing the alignment values to be 1 if we’re assuming that the panel is on a tracker. What’s more, we can measure the effect of tracking the azimuth only, the altitude only, or both at the same time!

We could add parameters to the code that would modify the behavior of the alignment function, but since I’m imagining that we’d run the code once with alignment enabled and once with a fixed panel, there’s not really much reason to make all of the code changes needed. Instead, I’ll introduce some global settings and modify just the alignment function:

TRACK_AZIMUTH = False
TRACK_ALTITUDE = False

def alignment(u, v):
    c = lambda a, b: np.cos(np.deg2rad(a - b))
    azimuth = 1 if TRACK_AZIMUTH else c(u.azimuth_deg, v.azimuth_deg)
    altitude = 1 if TRACK_ALTITUDE else c(u.altitude_deg, v.altitude_deg)
    return azimuth * altitude

In my case, enabling two-axis tracking improves the solar panel’s expected performance by a factor of about 2, which is enormous! Also, and I think quite interestingly, tracking the azimuth alone produces nearly all of the difference between a fixed panel and a two-axis tracking setup. I’d speculate that if a panel is to be located near the equator, tracking the sun’s altitude would provide most of the benefit of a tracking system. However, at latitudes further from the equator, tracking the azimuth probably matters more.

Full script by the end of orbit 5: panel-optimization-orbit5.py

Summary §

In this post I’ve presented a basic model for optimizing the output from a solar panel, and progressively refined the model to become more detailed. I hope that, along the way, my descriptions about the model have been edifying.

If you think any of this is noteworthy or find something that needs fixing, please let me know at hello@. Thanks!